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Abstract 

This  thesis  provides  the  groundwork  that  will  enable  development  of  a  lightweight, 
inexpensive,  aerodynamic,  and  broadband  antenna.  Whether  for  radar  or  communication, 
an  antenna  with  these  properties  would  be  a  force  multiplier  for  the  smaller,  limited  payload 
air  vehicles  the  United  States  Air  Force  will  pursue  in  the  coming  years. 

Several  microstrip  antennas  using  the  first  higher  order  mode  were  simulated  with 
the  Finite  Difference  Time  Domain  (FDTD)  method.  The  propagation  constant  of  each 
antenna  was  extracted  from  the  resulting  field  distribution  for  comparison  with  a  trans¬ 
verse  resonance  approximation,  measured  far-held  patterns,  and  other  simulated  antennas. 
Variations  of  the  geometry  were  explored  to  investigate  held  propagation,  improve  the  far- 
held  pattern,  and  improve  bandwidth.  A  simplified  fabrication  method  was  demonstrated 
that  shorten  production  time  and  improved  the  far-held  pattern. 
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FINITE  DIFFERENCE  TIME  DOMAIN  (FDTD)  ANALYSIS  OF  A 
LEAKY  TRAVELING  WAVE  MICROSTRIP  ANTENNA 


I.  Introduction 

The  antenna,  as  defined  in  The  Institute  of  Electrical  and  Electronic  Engineers  (IEEE) 
Standard  145-1983,  is  a  means  of  radiating  and  receiving  electromagnetic  waves.  In 
practice,  an  antenna  is  a  component  of  a  communication  system  that  provides  a  connec¬ 
tion  between  a  remotely  separated  sender  and  receiver.  Antennas  allow  the  sharing  of 
information  without  sharing  a  wire. 

The  need  for  more  information  sharing  between  individuals  is  making  our  world 
increasingly  wireless.  One  need  not  look  further  than  the  incredible  growth  of  cell  phones, 
pagers,  and  wireless  internet  connections.  Sometimes  the  sender  or  receiver  is  not  even 
a  person.  Unmanned  vehicles  that  operate  in  environments  too  expensive  or  too  hostile 
for  humans  can  be  remotely  controlled  using  an  antenna,  even  as  far  away  as  Neptune! 
Many  newer  automobiles  employ  antennas  to  see  behind  them  to  alert  the  driver  to  an 
obstruction. 

Aviation,  particularly  military  aviation,  uses  many  systems  whose  effectiveness  hinges 
upon  a  properly  functioning  antenna.  To  avert  mid-air  collisions,  aircraft  must  commu¬ 
nicate  in  real-time  with  ground  controllers  as  well  as  other  aircraft.  Landing,  which  is 
arguably  the  most  difficult  task  in  aviation,  is  made  relatively  simple  by  an  instrument 
landing  system  (ILS)  in  which  an  antenna  on  the  ground  talks  to  an  antenna  on  the  air¬ 
craft.  The  safety  and  enjoyment  of  flights  are  enhanced  by  weather  radar  systems  that 
allow  flights  to  be  diverted  around  storms.  Navigation  is  greatly  enhanced  from  systems 
that  use  antennas  to  communicate  with  satellites  or  scan  the  ground  with  radar.  Mission 
effectiveness  and  aircraft  survivability  are  increased  by  the  radar  warning  receivers,  which 
alert  a  pilot  to  the  invisible  dangers  of  enemy  aircraft  and  surface  to  air  missiles. 

Historically,  antennas  have  been  the  most  expensive  component  in  most  communica¬ 
tion  systems.  The  price  is  due  to  the  high  development  cost  of  designing  a  new  antenna. 
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The  electromagnetic  interactions  involved  in  the  operation  of  an  antenna  can  be  extremely 
complicated.  The  ability  to  analytically  characterize  many  new  antennas  is  prohibitively 
time  intensive.  Attention  is  paid,  instead,  to  numerically  approximating  the  structure, 
although,  an  accurate  model  often  requires  much  computing  power.  Adding  to  the  price 
tag,  an  antenna  must  be  manufactured  under  strict  tolerances  to  operate  as  designed. 

In  recent  years,  antenna  design  has  become  less  expensive.  Exponential  advances 
in  computing  resources  together  with  more  efficient  algorithms  have  created  a  climate 
where  elaborate  simulations  have  become  commonplace.  Improved  materials  and  new 
manufacturing  techniques  have  made  many  types  of  antennas,  such  as  microstrip,  easy  to 
fabricate. 

1.1  Problem  Statement 

The  U.  S.  Air  Force  has  a  need  for  a  conformal,  broadband  antenna  that  is  lightweight 
and  inexpensive.  A  conformal  antenna  is  more  aerodynamic  and  tends  to  have  a  lower 
radar  cross  section.  The  less  the  aircraft  weights,  the  more  payload  that  can  be  carried. 
Likewise,  the  less  the  aircraft  costs,  the  more  aircraft  that  can  be  purchased,  and  in  turn, 
the  more  payload  that  can  be  delivered.  A  broadband  antenna  reduces  the  need  for  multiple 
antennas,  which  further  reduces  the  cost  and  weight  of  the  aircraft.  Microstrip  inherently 
satisfies  all  of  these  requirements  except  bandwidth.  A  microstrip  antenna  operating  in  a 
leaky  traveling  wave  configuration  could  provide  the  bandwidth  needed. 

A  new  antenna  design  proposed  by  Dr.  Gary  Thiele  of  Analytic  Designs,  Inc.,  is 
seen  in  Figure  1.1.  The  design  is  based  on  work  by  Menzel  [23],  whose  antenna  can  be 
seen  in  Figure  1.2.  Menzel’s  antenna  uses  seven  slots  cut  from  the  conductor  along  the 
centerline  to  suppress  the  fundamental  mode  allowing  leaky  wave  radiation  via  the  first 
higher  order  mode.  Menzel’s  antenna  has  been  analyzed  by  a  host  of  researchers  over 
the  past  25  years  [2,9,10,21,24-26,28-30,34,40,42]  and  its  performance  is  known  and 
reproducible.  Instead  of  transverse  slots,  Thiele’s  antenna  uses  a  metal  wall  down  the 
centerline  of  the  antenna  to  block  the  fundamental  mode.  Symmetry  along  this  metal  wall 
invites  the  application  of  image  theory.  One  entire  side  of  the  antenna  is  now  an  image  of 
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Figure  1.1:  Thiele  Half  Width  (THW)  antenna. 


Matching  Stub 


Figure  1.2:  Menzel’s  original  antenna  [23]. 

the  other  side,  making  it  redundant  and  unneeded.  Thiele’s  resulting  antenna  is  half  the 
width  of  Menzel’s  antenna. 

This  thesis  modelled  Thiele’s  antenna  with  Finite  Difference  Time  Domain  (FDTD) 
techniques.  Different  geometries,  including  variations  of  curvature  of  the  conductor  strip, 
as  well  as  thickness  and  composition  of  the  substrate,  were  simulated.  Different  feeding 
methods  were  empirically  evaluated. 
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A  possible  application  of  this  new  antenna  could  include  an  inexpensive  replacement 
of  the  bulky  cavity-backed  spiral  antennas  used  for  the  great  majority  of  radar  warning 
receivers. 

1.2  Scope 

The  objective  of  this  research  was  to  characterize  the  propagation  constant  of  a  leaky 
mode  traveling  wave  antenna.  Understanding  how  the  propagation  constant  is  affected  by 
modifying  the  geometry  of  the  structure  is  vital  to  improving  the  far-held  pattern  and 
bandwidth  of  a  traveling  wave  antenna.  The  Finite  Difference  Time  Domain  (FDTD) 
method,  which  is  a  computational  electromagnetic  (CEM)  technique,  was  used  to  simulate 
different  antenna  geometries.  The  propagation  constant  was  extracted  from  the  held  dis¬ 
tribution  resulting  from  the  FDTD  simulation.  A  design  was  fabricated  and  tested  that 
incorporated  improvements  from  simulations. 

1.3  Resources 

Code  was  written  in  Matlab  starting  from  a  program  written  by  Keely  Willis  and 
Dr.  Susan  Hagness  of  the  University  of  Wisconsin  Computational  Electromagnetics  Labo¬ 
ratory.  Fabrication  required  computer-aided  design  (CAD)  software,  microstrip  material, 
coaxial  feed  materials,  and  milling  equipment.  Measurements  required  a  network  analyzer, 
an  antenna  test  range,  and  associated  supplies. 

1-4  Overview 

This  thesis  presents  theories  related  to  leaky  wave  microstrip  antennas,  methods  of 
modelling  such  structures  with  Finite  Difference  Time  Domain  (FDTD),  and  both  ex¬ 
perimental  and  simulated  results  of  several  antennas  based  on  Thiele’s  proposed  design. 
Chapter  II  provides  a  literature  review  of  work  related  to  leaky  wave  microstrip  antennas. 
Chapter  III  is  a  discussion  of  the  FDTD  method.  Chapter  IV  describes  development  of 
the  FDTD  simulation.  Chapter  V  presents  results  and  analysis  of  tested  antenna  designs. 
Finally,  Chapter  VI  states  conclusions  and  suggests  recommendations  for  future  research. 
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II.  Background  on  Microstrip  Leaky  Traveling  Wave  Antennas 
2. 1  Traveling  Waves 

Much  like  the  waves  produced  by  dropping  a  stone  into  a  pond,  a  traveling  wave 
is  an  electromagnetic  disturbance  that  propagates  with  a  constant  phase.  Since  the  wave 
maintains  a  constant  phase,  the  wave  travels  with  a  constant  phase  velocity,  vp.  Obeying 
Maxwell’s  equations  dictates  that  this  wave  must  also  satisfy  the  Helmholtz  wave  equation. 
For  instance,  the  source-free  vector  wave  equation  of  a  magnetic  field  in  a  lossless  medium 
is: 


(V2  +  k2)il  =  0  (2.1) 

The  solution  to  Equation  (2.1)  for  a  complex  magnetic  field  plane  wave  polarized  in  the 
y  direction  and  traveling  in  the  +x  direction  is: 

H  =  yHyej{u)t~kx)  (2.2) 

u j  is  the  angular  frequency  of  the  excitation  source.  The  propagation  constant,  k,  is  com¬ 
posed  of  a  real  portion,  (3 ,  which  is  called  the  phase  constant,  and  an  imaginary  component, 
a,  called  the  attenuation  constant: 


k  =  f3  —  ja 


(2.3) 


The  phase  velocity  is  given  by: 


oj  1 

~P  V^i1 


(2.4) 


where  e  and  g  are  the  constitutive  parameters  of  the  material  in  which  the  wave  is  propa¬ 
gating. 
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2.2  Traveling  Wave  Antennas 

Most  antennas  operate  in  a  resonance  configuration  by  supporting  standing  waves  of 
the  currents,  voltages,  and  fields  along  their  length,  which  is  analogous  to  the  vibrations  of 
a  violin  string.  Standing  waves  can  be  described  as  a  superposition  of  two  traveling  waves 
propagating  in  opposite  directions.  Dipoles  and  microstrip  patches  are  common  examples 
of  antennas  in  which  two  waves  of  equal  amplitude  that  are  180°  out  of  phase  travel  in 
opposite  directions  along  their  length.  The  current  distribution  is  a  standing  wave  with 
nulls  at  the  ends,  while  the  voltage  distribution  is  a  standing  wave  with  maxima  at  the 
ends.  The  standing  waves  result  when  energy  is  excited  at  one  end,  travels  the  length 
of  the  antenna,  reflects  off  of  the  opposite  end,  and  returns  toward  the  feed.  A  result  of 
standing  waves  is  limited  bandwidth  due  to  the  destructive  interference  of  the  waves  [41]. 

A  solution  to  this  problem  is  an  antenna  that  uses  a  single  traveling  wave.  The 
standing  wave  is  prohibited  by  eliminating  the  wave  in  the  return  direction,  which  can  be 
accomplished  by  placing  an  impedance  at  the  far  end  to  dissipate  any  remaining  energy. 
Alternatively,  the  antenna  can  be  made  long  enough  so  as  to  radiate  nearly  all  of  the 
energy  before  the  forward  wave  from  the  feed  reaches  the  opposite  end,  however,  this 
requires  the  antenna  to  become  unreasonably  long.  A  general  rule  is  to  design  the  length, 
L,  to  radiate  90%  of  the  applied  power,  P ,  and  dissipate  the  remaining  10%  with  a  load, 
as  in  Equation  (2.5): 


P(L) 

m 


=  e~2aL 


=  0.1 


(2.5) 


Equation  (2.5)  directly  leads  to  Equation  (2.6),  which  relates  normalized  a  to  normalized 
length. 


L  -0.183 

Ao  is  the  excitation  wavelength  in  free  space  and  ko  is  the  wavenumber  in  free  space. 

The  simplest  example  of  a  traveling  wave  antenna  is  a  long  wire  supported  a  distance, 
h,  above  the  ground,  also  known  as  a  Beverage  antenna.  Figure  2.1  shows  a  Beverage 


6 


Figure  2.1:  A  long  wire  antenna,  or  Beverage  antenna,  illustrating  the  effects  of  a 

reflected  traveling  wave  [4]. 

antenna  excited  by  a  forward  current,  If.  Without  a  matched  resistor,  Rl,  a  backward 
current,  R,  is  created  by  the  reflection  of  If  at  the  far  end.  A  forward  lobe  is  produced 
by  If  and  a  backward  lobe  is  produced  by  R.  In  Figure  2.1,  when  the  impedance  of  Rj_, 
is  matched  to  the  end  of  the  wire  antenna,  R  is  eliminated,  which  removes  the  backward 
lobe  associated  with  the  returning  traveling  wave.  Since  an  antenna  that  is  designed  to 
operate  with  a  traveling  wave  will  have  uniform  patterns  of  current  and  voltage,  it  behaves 
like  a  resistive  circuit  component.  Matching  the  feed  line  to  the  antenna  is  simple  since 
the  antenna  has  a  purely  real  input  impedance. 

2.3  Antenna  Characteristics 

The  current  can  be  approximated  as  a  finite  line  source  that  radiates  as  an  outward 
propagating  spherical  wave.  The  spherical  wave  can  be  approximated  as  a  plane  wave  in 
the  far-held.  Equation  (2.7)  determines  the  far-held  source  (or  observation)  point  distance, 
r,  at  which  the  phase  varies  by  no  more  than  over  half  of  the  length,  L,  as  depicted  in 
Figure  2.2. 
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Spherical 


Figure  2.2:  The  far- field  begins  at  distance  r  from  the  antenna  at  which  the  spherical 

wavefront  varies  by  less  than  over  half  of  the  antenna’s  longest  dimension. 
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(2.7) 


The  current  distribution  along  the  length  of  the  traveling  wave  antenna  is: 


/  =  Ioe~jkz 


(2.8) 


The  vector  potential,  A(x,y,z),  of  an  antenna  with  /  is  then: 

„  fi  fL  I(z')e-ik°r'  ,  , 

A  (x,y,z)  =  zAz  =  zj-  - - - dz 

4tt  J o  r' 


(2.9) 
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After  making  the  far-held  approximation,  Ar  =  0,  the  spherical  coordinate  equivalent 
vector  potential  is: 


A (r,  0,  0)  =  0Ae  =  -§AZ  sm6  =  _^/pe  ^°rsm^  fL  ej(kocos0-k)z' dz>  (2.10) 

4tt  r  J  o 


The  electric  held,  E,  due  to  A  follows: 


E=-jnA  =  fLeKM,co„-kWd2, 

4vrr  J0 


(2.11) 


Solving  the  integral  gives: 


= 


jcc/j/oe  ifc°r  sin  0 
47T?’ 


gj(ko  cos  0—k)L  ^ 

j  (&0  cos  6  —  k) 


(2.12) 


which  can  be  separated  into  an  element  factor,  Eo,  and  a  pattern  factor  as  in  Equa¬ 
tion  (2.13): 


Eg  =  Eq- 


sin  O 


n 


(2.13) 


where  the  element  factor  is: 


p  jufiloL  sin  0e  jk°r  jk(knCOSe-k) 

&0  =  - ~A -  Z  2  V  ' 

47 rr 


(2.14) 


and  the  pattern  factor  argument  is: 


12  =  —  (kocosO  —  k)  (2.15) 

The  observation  point  is  assumed  in  the  x  —  z  plane,  where  0  =  0.  The  pattern  factor  is 
a  maximum  when  its  argument,  fl,  equals  zero.  This  forms  the  main  lobe  at  6m,  given  by 
Equation  (2.16). 
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6.7  GHz  a  /k  =  0.036  P  /k  =  0.67 


Figure  2.3:  The  far- field  electric  field  pattern  of  a  traveling  wave  antenna  with  aperture 
three  wavelengths  long. 


17  =  0  =  k$  cos  9m  —  k 

9m  =  cos-1  j-  (2.16) 

«o 

The  magnetic  field,  H.  is  found  by: 

H  =  —  r  x  E  (2.17) 

Vo 

where  t/q  is  the  impedance  of  free  space. 

Equations  (2.13)-(2.15)  produce  the  pattern  seen  in  Figure  2.3  for  the  Thiele  Half 
Width  (THW)  antenna  of  the  dimensions  in  Figure  1.1,  which  is  roughly  three  Xg  long. 
Xp  is  the  wavelength  of  the  traveling  wave.  Notice  the  main  beam  is  not  directed  at 
the  angle  predicted  by  Equation  (2.16)  of  44.9°.  The  length  must  be  greater  than  five 
wavelengths  long  and  9m  >  20°  for  Equation  (2.16)  to  provide  an  approximation  within 
1°  [41].  Figure  2.4  is  a  pattern  of  the  same  antenna  lengthened  to  approximately  5  Xg  long. 


10 


6.7  GHz  a  /k  =  0.036  P  /k  =  0.67 


Figure  2.4:  The  far-field  electric  field  pattern  of  the  same  antenna  in  Figure  2.3  with 

aperture  increased  to  five  wavelengths  long. 

The  far-field  pattern  shown  in  Figures  2.3  and  2.4  only  includes  the  forward  traveling 
wave.  Of  course,  a  finite  length  antenna  will  have  a  finite  amount  of  energy  reflecting  from 
the  end,  returning  in  the  opposite  direction,  and  producing  a  backward  lobe.  Figure  2.5 
shows  the  effect  of  summing  four  traveling  waves,  using  reflection  coefficients  of  -1  at  both 
ends. 

When  transmitting,  the  structure  of  the  antenna  is  constructed  to  radiate  as  much 
energy  as  possible  into  free  space  for  a  given  applied  current.  Gain  is  a  typical  measure 
of  the  efficiency  of  an  antenna  to  radiate  in  a  desired  direction.  In  Equation  (2.18),  gain, 
G,  is  the  dimensionless  ratio  of  the  radiation  intensity  at  a  given  direction,  U(9,  (j>),  to  the 
radiation  intensity  obtained  if  the  input  power,  P-m ,  were  radiated  isotropically. 


G 


4irU  ( 9 ,  4>) 


Pit 


(2.18) 


There  are  two  basic  means  of  describing  bandwidth,  which  is  the  range  of  frequencies 
of  acceptable  performance.  The  bandwidth  of  broadband  antennas  is  often  given  by  the 
ratio  of  the  highest  frequency,  fn,  to  the  lowest  frequency,  fj,.  For  example,  an  antenna 
that  operates  on  the  band  2-20  GHz  would  have  a  10:1  bandwidth.  Narrowband  antennas 
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6.7  GHz  a  /kQ=  0.036  p  /kQ=  0.67 


Figure  2.5:  The  far- field  pattern  generated  by  a  traveling  wave  with  three  reflections  at 

6.7  GHz. 

usually  express  their  bandwidth  as  a  percentage  of  the  center  frequency,  fc.  An  antenna 
that  operates  in  the  band  7-9  GHz  would  have  the  following  bandwidth: 


Bandwidth  =  —  ^ L  =  — — —  =  25%  (2.19) 

Jc  8 

A  result  of  the  Lorentz  reciprocity  theorem,  which  is  derived  from  Maxwell’s  equa¬ 
tions,  is  that  the  gain  pattern  of  an  antenna  is  identical  whether  it  is  transmitting  or 
receiving,  as  long  as  it  is  surrounded  by  a  linear,  isotropic  medium.  Simulation  of  trans¬ 
mission  alone  is  all  that  is  needed  to  fully  characterize  an  antenna. 


2-4  Propagation  Modes 

A  mode  is  a  particular  configuration  of  the  fields  in  a  transmission  line.  A  Transverse 
Electric  and  Magnetic  (TEM)  mode  is  a  field  distribution  at  which  both  the  electric,  E, 
and  magnetic,  H,  field  intensities  are  contained  in  an  equiphase  plane  that  is  independent 
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of  time  [3].  This  can  also  be  described  as  both  E  and  H  being  transverse  to  the  direction 
of  propagation.  A  higher  order  mode  is  a  mode  in  which  either  E  or  H  has  a  component 
in  the  direction  of  propagation.  For  example,  a  TEn,  or  Transverse  Electric  mode,  wave 
traveling  in  the  z  direction  has  a  vector  component  of  H  in  the  z  direction  while  all  vector 
components  of  E  are  transverse  to  z.  The  n  signifies  that  there  are  an  infinite  number  of 
these  modes  that  can  be  found  using  Maxwell’s  equations.  Likewise,  TM„  signifies  the  nth 
higher  order  Transverse  Magnetic  mode.  While  many  modes  may  be  possible  in  a  structure, 
the  majority  of  the  energy  tends  to  dominate  the  lowest  mode  that  can  propagate. 

In  a  given  transmission  line  at  a  given  excitation  frequency,  the  modes  that  can  exist 
are  determined  by  each  mode’s  cutoff  frequency.  A  mode’s  cutoff  frequency  depends  upon 
the  dimensions  of  the  structure  and  the  medium  inside  the  structure  [3] .  The  modes  which 
have  a  cutoff  frequency  equal  to  or  smaller  than  the  excitation  frequency  are  supported 
and  all  others  will  quickly  decay. 

2. 5  Microstrip 

2.5.1  Physical  Characteristics.  Microstrip  was  first  proposed  by  Deschamps  in 
1953  at  the  3rd  Air  Force  Symposium  on  Antennas  [18],  and  it  remained  an  academic 
novelty  for  nearly  20  years.  Advances  in  materials  and  manufacturing  processes  in  the 
1970’s  made  its  production  feasible.  Microstrip  antenna  technology  has  been  the  most 
rapidly  developing  topic  in  antennas  during  the  last  twenty  years  [32],  Microstrip  is  an 
open  structure  that  consists  of  a  very  thin  metallic  strip  of  a  width,  w,  separated  from 
a  ground  plate  by  a  dielectric  sheet  called  substrate  (Figure  2.6).  The  thickness  of  the 
conductor,  t,  is  much  less  than  a  wavelength.  The  height  of  the  substrate,  h,  is  usually 
very  thin  compared  to  the  wavelength  (.0003A  <  h  <  0.05A)  [31].  The  substrate  is  designed 
to  have  a  known  relative  permittivity,  er ,  that  is  homogeneous  within  specified  temperature 
limits. 

Figure  2.7  illustrates  some  of  the  many  methods  to  feed  a  microstrip  antenna.  The 
antenna  can  be  excited  directly  by  a  microstrip  line,  by  a  coaxial  cable,  or  a  combination  of 
the  two.  The  antenna  can  also  be  fed  from  a  microstrip  line  without  direct  contact  through 
electromagnetic  coupling.  Feeding  by  electromagnetic  coupling  through  an  aperture  in  the 
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Figure  2.6:  Geometry  of  a  microstrip  transmission  line. 

ground  plane  tends  to  improve  bandwidth.  To  maximize  efficiency,  the  impedance  of  the 
feed  must  be  matched  to  the  input  impedance  of  the  antenna.  There  are  a  variety  of  stubs, 
shunts,  and  other  devices  used  for  matching. 

2.5.2  Advantages  and  Disadvantages.  The  advantages  of  microstrip  lie  in  its 
physical  characteristics.  Circuits  made  with  microstrip,  especially  antennas,  can  be  essen¬ 
tial  to  aircraft,  spacecraft,  satellite,  and  missile  applications.  Not  only  do  these  circuits 
have  an  inherently  aerodynamic  profile,  but  they  can  also  be  conformable  to  a  surface.  Mi¬ 
crostrip  is  lighter  and  smaller  than  parabolic  dishes  and  waveguide  arrays,  and  generally 
cheaper  and  easier  to  manufacture  and  install  [4, 18] .  Mass  production  printed  circuit  and 
chemical  etching  technology  have  led  to  very  low  fabrication  costs  [18]. 

Compared  to  other  microwave  antennas,  the  major  disadvantages  of  microstrip  are 
lower  gain  and  very  narrow  bandwidth  [18].  Microstrip  characteristically  has  low  efficiency 
and  low  power  handling  ability  [4,32].  In  addition,  antennas  made  with  microstrip  typically 
have  poor  polarization  purity  and  poor  scan  performance  [18]. 

2.5.3  Hybrid  Modes.  Operating  above  the  cutoff  frequency,  the  field  lines  of 
microstrip  extend  throughout  the  substrate  as  well  as  into  the  free  space  region  above  the 
substrate,  as  seen  in  Figure  2.8.  The  phase  velocity  of  the  field  in  the  free  space  surrounding 
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Figure  2.7:  Exploded  view  of  microstrip  feed  techniques  for  a  square  patch:  (a)  Direct, 
(b)  Electromagnetic  coupling,  (c)  Coupling  through  an  aperture  in  the  ground  plane,  (d) 
Direct  coaxial  connection  [18]. 


the  structure  is  the  speed  of  light,  c,  and  the  phase  velocity  of  the  field  in  the  substrate  is 
given  by  Equation  (2.20). 


vn  = 


(2.20) 


This  difference  in  phase  velocity  at  the  interface  between  the  substrate  and  free  space  makes 
the  TEM  mode  impossible.  Instead,  the  fundamental  mode  for  microstrip  is  a  quasi-TEM 
mode,  usually  annotated,  EHo,  in  which  both  the  electric  and  magnetic  fields  have  a 
component  in  the  direction  of  propagation.  Likewise,  a  higher  order  mode  in  microstrip  is 
not  purely  TE  or  TM,  but  a  hybrid  combination  of  the  two.  The  nth  higher  order  mode 
is  termed  the  EHn  mode. 

The  fundamental  mode  of  microstrip,  as  seen  in  Figure  2.8,  does  not  radiate  since 
the  fields  produced  do  not  decouple  from  the  structure.  If  the  fundamental  mode  is  not 
allowed  to  propagate,  the  next  higher  order  mode  will  dominate.  Figure  2.9  shows  the 
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Figure  2.8:  Field  pattern  associated  with  the  fundamental  mode  of  microstrip,  EHq. 


►  E-field  lines 


Figure  2.9:  Field  pattern  associated  with  the  first  higher  order  mode  of  microstrip,  EHi. 

fields  due  to  the  first  higher  order  mode,  EHi.  A  phase  reversal,  or  null,  appears  along  the 
centerline,  allowing  the  fields  to  decouple  and  radiate. 

2.6  Propagation  Mechanisms 

Higher  order  modes  on  microstrip  transmission  lines  can  be  described  as  exhibiting 
three  distinct  propagation  mechanisms  above  the  cutoff  frequency,  fc:  bound  wave,  surface 
wave,  and  leaky  wave.  These  mechanisms  can  be  described  by  analyzing  the  general 
dispersion  relation  of  Equation  (2.21): 
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Figure  2.10:  The  typical  normalized  attenuation  constant,  ,  and  phase  constant,  in 
the  direction  of  propagation  of  the  first  higher  order  mode,  EHi.  There  are  four  frequency 
regions  associated  with  propagation  regimes:  Reactive,  Leaky,  Surface,  and  Bound. 


k  —  kx  T  ky  T  kz 


(2.21) 


where  the  complex  wavenumbers  are  composed  of  a  phase  constant,  (3,  and  an  attenuation 
constant,  a: 


kx  —  flx  JCXx 

ky  =  fly  —  jay 

kz  =  flz  -  jaz  (2.22) 

Following  the  orientation  of  Figure  2.6,  kx  is  the  wavenumber  in  the  direction  of 
propagation.  Figure  2.10  is  a  typical  plot  of  kx  over  each  of  the  propagation  mechanism 
regions.  Below  the  cutoff  frequency,  ax  is  quite  large  causing  the  transmission  line  to 
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function  as  a  reactive  load.  Unsurprisingly,  this  region  is  frequently  called  the  reactive 
region. 

As  the  frequency  is  raised  to  fc,  /3X  =  ax  resulting  in  the  commonly  used  cutoff 
frequency  expression:  kx  =  0.  The  dispersion  relation  at  this  point  is: 


k 2 


Py  +  Pz 

L02£fJ, 


(2.23) 


which  can  be  manipulated  to  find  the  cutoff  frequency,  fc: 


fc  = 


(2.24) 


Directly  above  cutoff  is  the  leaky  region  where  energy  begins  to  propagate  down 
the  transmission  line  as  (3X  grows  larger  than  ax.  Field  losses  due  to  ax  are  not  ohmic. 
Although  surface  waves  are  present,  losses  in  the  leaky  wave  region  are  mostly  due  to 
the  radiation  of  energy  leaking  from  the  microstrip  [2],  A  purely  real  kz  allows  power  to 
radiate,  or  couple  per  unit  length  from  the  structure  into  free  space  [13]. 

The  leaking  wave  travels  away  from  the  antenna  at  an  angle,  9,  measured  from  the 
endfire  direction  as  seen  in  Figure  2.11.  At  the  lower  limit  of  the  leaky  transition  region, 
the  beam  is  nearly  broadside.  0  decreases  as  the  frequency  increases.  Simple  geometry 
shows  the  relation  between  frequency  and  8  in  Equation  (2.25),  which  is  equivalent  to 
Equation  (2.16).  The  main  beam  approaches  endfire  as  the  frequency  approaches  the 
leaky  region  upper  limit,  at  which  point  px  =  ko- 


-1 

'  Px " 

-1 

PxC 

cos 

-ko- 

=  cos 

.  tv  . 

(2.25) 


For  any  frequency  in  the  leaky  region,  there  is  a  curious  property  whereby  the  fields 
due  to  the  leaky  wave  actually  increase  moving  away  from  the  structure.  The  reason 
for  this  is  simply  that  less  energy  leaks  per  unit  length  as  the  wave  travels  down  the 
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broadside 


Figure  2.11:  The  direction  of  the  main  beam  radiation  occurs  at  an  angle  created  by  (3X 
and  ko. 


broadside 


endfire 


Figure  2.12:  The  intensity  of  fields  radiating  from  a  leaky  wave  antenna  increase  moving 
away  from  the  structure. 


structure.  In  Figure  2.12,  the  strength  of  the  leaking  field  is  depicted  by  the  thickness  of 
the  line.  The  fields  increase  exponentially  to  a  distance  above  the  antenna,  zmax ,  given  by 
Equation  (2.26),  and  then  quickly  decay  [29]. 


Zmax(x)  =  x  tan  9  (2.26) 

where  x  is  the  distance  from  the  source  feed  in  the  x  direction.  Rotating  angle  9  about 
the  x  axis  of  the  antenna  forms  a  bowl-shaped  main  lobe. 
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Leaky  waves  are  considered  fast  waves  since  the  phase  front  travels  faster  than  the 
speed  of  light.  A  leaky  wave  is  generally  an  efficient  means  to  radiate  energy,  but  it  is 
typically  narrow-banded  since  it  occurs  only  for  frequencies  in  a  transition  region,  as  seen 
in  Figure  2.10. 

Above  the  frequency  at  which  fix  =  ko,  az  increases  causing  the  leaky  wave  to 
attenuate.  The  surface  wave,  however,  remains.  Surface  waves  are  unbounded  losses 
that  are  supported  by  a  dielectric  layer  on  a  ground  plate.  The  magnitude  of  a  surface 
wave  decays  exponentially  as  it  travels  away  from  the  transmission  line  in  the  —z  direction 
since  kz  has  a  large  imaginary  component,  therefore,  these  waves  are  associated  only  with 
the  surface  of  a  structure.  Surface  waves  can  emanate  from  discontinuities  in  the  physical 
structure  that  disrupt  and  decouple  the  bound  waves  from  the  surface.  Discontinuities 
can  include  ends,  corners,  feeding  structures,  and  also  curvature.  Generally  these  surface 
waves  are  classified  as  slow  waves  since  their  phase  velocity  is  less  than  the  speed  of  light. 
Surface  waves  are  the  means  by  which  desired  coupling  takes  place  in  certain  microwave 
circuits  such  as  coupled  line  filters,  however,  more  often,  they  tend  to  complicate  both 
antenna  and  circuit  design  by  introducing  additional  loss  and  unwanted  coupling  [6]. 

Just  like  leaky  waves,  surface  waves  travel  outward  from  the  microstrip  at  an  angle,  9. 
Referring  to  Figure  2.13,  the  wavenumber  of  the  surface  wave,  ks,  has  a  component  in  the 
x  direction,  (3X,  and  a  component  in  the  y  direction,  (3y.  The  resulting  dispersion  relation, 
Equation  (2.27),  shows  that  j3x  <  ko^/ef  for  surface  waves  to  exist: 


Py  =  k2s-Pl 

(3y  =  yJkfer-(%  (2.27) 

Above  the  frequency  at  which  /3X  =  ks,  ay  increases  causing  the  surface  wave  to 
attenuate.  The  field  waves  are  then  bound,  to  the  immediate  vicinity  of  the  structure.  This 
regime  is  typically  the  state  at  which  most  microwave  circuits  are  designed  to  operate.  The 
operating  frequency  is  generally  chosen  to  be  20%  above  fc  to  ensure  leaky  and  surface 
wave  losses  are  not  present. 
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Figure  2.13:  The  surface  wavenumber,  ks,  has  an  x  component,  f3x,  and  a  y  component, 

Py 

2.7  Menzel’s  Original  Antenna 

Microstrip  structures  do  not  radiate  for  the  fundamental  mode,  therefore,  a  higher 
order  mode  must  be  excited.  In  1979,  Menzel,  published  the  first  account  of  a  traveling 
wave  microstrip  antenna  that  used  a  higher  order  mode  to  produce  leaky  waves  [23] .  This 
method  of  producing  radiation  by  exciting  higher  order  modes  in  a  transmission  line  has 
been  documented  since  the  1950’s  [11].  By  the  1970’s,  rectangular  waveguides,  circular 
waveguides,  and  coaxial  cables  were  in  use  as  leaky  traveling  wave  antennas.  However, 
until  Menzel,  the  jump  to  microstrip  had  not  been  made. 

By  looking  at  a  cross  section  of  microstrip  excited  in  the  fundamental  mode,  the 
E  field  is  strongest  in  the  center  and  tapers  off  to  zero  at  the  sides,  as  depicted  in  Fig¬ 
ure  2.8.  If  the  electric  field  down  the  centerline  is  suppressed,  the  fundamental  mode  will 
be  prohibited,  forcing  the  energy  to  propagate  at  the  next  higher  mode,  EHi.  As  seen  in 
Figure  2.9,  EHi  mode  causes  E  to  be  strongest  at  the  edges.  Menzel  attempted  to  force 
the  EHi  mode  using  several  means.  Feeding  two  equal  magnitude  waves  180°  out  of  phase 
with  a  “T”  or  “Y”  feed  produced  EHi  as  desired,  but  did  not  fully  eliminate  the  EHo 
mode.  Metal  posts,  known  as  vias,  inserted  down  the  centerline  eliminated  the  fundamen- 
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tal  mode  and  produced  the  anticipated  radiation.  A  drawback  of  using  vias  was  difficulty 
constructing  an  antenna  that  was  no  longer  planer.  Easier  to  produce  and  providing  an 
even  better  response  was  given  using  transverse  slots  down  the  centerline  (Figure  1.2).  The 
multiple  feeds  were  not  necessary  to  produce  the  EHi  mode  when  the  fundamental  mode 
was  suppressed. 

Menzel  demonstrated  that  the  beam  angle  can  be  predictively  steered  by  input  fre¬ 
quency  if  the  electrical  length  of  the  antenna  is  at  least  3A.  If  the  length  is  less  than  3A, 
too  little  of  the  incident  wave  is  being  radiated  and  a  resonance  standing  wave  pattern  is 
forcing  the  beam  toward  broadside. 

Qualitative  analysis  shows  that  the  beamwidth  of  Menzel’s  antenna  is  not  frequency 
dependent,  however,  it  is  inversely  related  to  length.  The  3  dB  beamwidth  approaches  10° 
for  electrical  length  of  over  6A  and  approaches  nearly  90°  for  fractions  of  a  wavelength. 

Menzel’s  gain  varied  from  7  dB  for  l  =  0.2A  to  14  dB  for  l  =  4A.  7  dB  is  comparable 
to  a  similar  sized  resonant  antenna.  An  antenna  longer  than  l  =  4A  would  have  an  even 
higher  gain  as  the  radiation  aperture  increases. 

2.8  Analysis  of  Menzel’s  Work 

Lee  notes  that  Menzel  assumed  that  his  antenna  should  radiate  simply  because  the 
phase  constant  due  to  his  operating  frequency  was  less  than  ko  [20].  If  Menzel  had  con¬ 
sidered  the  complex  propagation  constant,  he  would  have  realized  that  his  antenna  was 
operating  in  a  leaky  regime.  The  length  would  need  to  be  roughly  220  mm,  or  more  than 
twice  as  long  as  his  design,  to  radiate  at  90%  efficiency.  Radiation  patterns  in  Menzel’s 
paper  clearly  show  the  presence  of  a  large  backlobe  due  to  the  reflected  traveling  wave. 

Errnert,  who  conducted  a  mode  matching  analysis  of  a  similar  microstrip  structure, 
did  not  investigate  normalized  phase  constants  below  j3  =  ko.  Errnert  specifically  refused 
to  include  leaky  modes  citing  Marcuvitz,  who  declared  these  waves  not  part  of  the  complete 
set  of  eigenmodes  [10,22].  Marcuvitz  rejected  these  modes  since  their  magnitude  appears 
to  increase  indefinitely  away  from  the  antenna  and,  hence,  do  not  satisfy  the  radiation 
conditions  at  infinity.  As  mentioned  earlier,  Figure  2.12  shows  that  the  field  strength 
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Figure  2.14: 


Predictions  of  normalized  phase  constant  for  Menzel’s  antenna  [29]. 


increases  exponentially,  but  it’s  maximum  is  at  a  finite  distance  given  by  Equation  (2.26). 
Bagby,  Grimm,  Nyquist  et  al  showed  that  while  leaky  modes  are  indeed  nonspectral,  a 
Steepest  Descent  Contour  (SDC)  integration  can  be  used  to  demonstrate  that  these  modes 
approximate  the  continuous  radiation  spectrum  [2,12], 

Oliner  used  a  transverse  resonance  formulation  to  also  show  that  the  leaky  region 
should  not  be  neglected.  The  normalized  phase  constant  in  Figure  2.14  predicted  by 
Menzel’s  calculations,  which  ignores  the  leaky  region,  is  shown  atop  that  predicted  by 
Oliner,  which  accounts  for  the  leaky  regime  [23,29].  Oliner  demonstrated  the  reason  for 
the  less  than  expected  bandwidth  experienced  by  Menzel  was  due  to  his  neglecting  of  the 
leaky  region.  Oliner’s  Steepest  Descent  Plane  analysis  predicts  that  the  useful  bandwidth 
of  the  higher  order  modes  can  be  quite  substantial. 

Michalski  and  Zheng  devised  a  rigorous  Mixed  Potential  Electric  Field  Integral  Equa¬ 
tion  (MPEFIE)  solution  that  was  applicable  to  an  arbitrary  cross  section.  For  Men¬ 
zel’s  antenna,  Michalski  and  Zheng’s  MPEFIE  dispersion  curves  are  in  agreement  with 
Oliner’s  [26,29]. 
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The  four  analytical  methods  mentioned,  Steepest  Descent,  Transverse  Resonance, 
Mode  Matching,  and  MPEFIE,  all  provide  much  needed  insight  into  the  physical  nature 
of  the  leaky  wave  phenomenon,  however,  they  all  have  the  very  serious  drawback  that  the 
solution  must  be  reformulated  for  each  geometry.  The  formulation  of  these  methods  is 
quite  complicated  and  time-intensive. 

Sheen  et  al  introduced  an  S-parameter  extraction  technique  to  determine  the  com¬ 
plex  propagation  constant  through  network  analyzer  measurements.  The  method  has  been 
validated  with  the  spectral  domain  approach  of  [2, 12].  Near-field  probing  measurements 
conducted  by  Thiele  show  promise  to  directly  measure  kx.  The  problem  with  any  mea¬ 
surement  is  that  the  fields  must  necessarily  be  disturbed  by  the  testing  equipment. 

In  lieu  of  analytical  methods  and  measurement  methods,  this  work  attempted  to 
develop  a  numerical  simulation  of  a  leaky  microstrip  antenna  that  would  be  easy  to  modify 
for  new  geometries.  The  Finite  Difference  Time  Domain  approach  was  chosen  since  the 
technique  directly  solves  for  the  fields  in  the  time  domain  using  Maxwell’s  equations. 
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III.  Finite  Difference  Time  Domain 


Finite  Difference  Time  Domain  (FDTD)  is  a  computational  electromagnetic  (CEM) 
technique  that  directly  solves  the  differential  form  of  Maxwell’s  equations, the  curl 
equations,  in  the  time  domain  using  a  discretized  space-time  grid.  Compared  to  an  integral 
equation  solution  of  Maxwell’s  equations,  such  as  Method  of  Moments  (MoM),  FDTD  offers 
the  benefits  of  no  linear  algebra,  well-understood  error  sources,  impulse  behavior  that  is 
treated  naturally,  nonlinear  behavior  that  is  treated  naturally,  enhanced  visualization  of 
the  wave  interactions,  and  a  systematic  approach  that  does  not  require  reformulations  of 
integral  equations  for  each  new  structure. 

Yee  first  proposed  using  finite  differencing  for  electromagnetic  problems  in  1966  [43]. 
Taflove  conducted  much  of  the  initial  research  into  Yee’s  method  in  the  1970’s  and  coined 
the  term  FDTD  in  1981  [39].  Although  there  are  other  methods  of  gridding  the  computa¬ 
tional  space,  the  grid  in  this  project  follows  the  Yee  Cell,  as  seen  in  Figure  3.1.  The  Yee 
cell  uses  rectangular  coordinates  and  places  each  of  the  components  of  the  magnetic  field 
a  half  grid  space  apart  from  the  orthogonally  directed  electric  field  components.  The  same 
must  necessarily  be  true  for  the  electric  field  with  respect  to  the  magnetic  field.  Not  only 
are  the  E  and  H  fields  a  half  of  a  cell  width  apart,  but  they  are  also  updated  a  half  of 
a  time  step  apart.  The  Yee  method  involves  a  leap  frog  approach  in  which  the  E  Held  is 
updated,  then  the  H  field,  then  E,  then  H,  and  so  on.  One  time  step  is  counted  after  both 
E  and  H  have  been  updated  once. 

3. 1  Formulation 

Yee’s  Finite  Difference  scheme  uses  central  differencing  that  is  second  order  accurate 
in  both  space  and  time.  Starting  with  Maxwell’s  curl  equations: 

—  =  -VxE-M  (3.1) 


<9D 

~dt 


V  x  H- J 


(3.2) 
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Figure  3.1:  The  Yee  Cell  [39,43]. 


Noting  that: 


J 


=  J 


source 


+  ctE 


(3.3) 


M  =  Msource  +  a*  H 


(3.4) 


Equations  (3.1)  and  (3.2)  can  be  manipulated  into  a  more  usable  form  that  is  applicable 
for  all  linear,  isotropic,  nondispersive,  lossy  materials: 


<9H  1  1 

— —  = - V  X  E - (Msource  +  CT*H)  (3-5) 

at  n  fx 

<9E  1  1 

— — -  =  -V  X  H - (j source  +  °"E)  (3-6) 

at  e  e 

Considering  only  the  source-free  case  and  separating  Equations  (3.5)  and  (3.6)  into  the 
vector  components  yields: 
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(3.7) 


dHx 

1 

\dEy 

dEz 

dt 

tL 

dz 

dy 

dHy 

1 

\dEz 

dEx 

dt 

h 

dx 

dz 

dHz 

1 

dEx 

dEy 

dt 

M 

dy 

dx 

dEx 

1 

\dHz 

dHy 

dt 

e 

dy 

dz 

dEy 

1 

\dHx 

dHz 

dt 

£ 

dz 

dx 

9EZ 

1 

\  dHy 

dHx 

dt 

£ 

dx 

dy 

The  finite  difference  implementation  of  Equation  (3.7)  is  as  follows: 


(3.8) 

(3.9) 

(3.10) 

(3.11) 

(3.12) 
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(3.13) 


The  convention  used  here  is  to  update  the  E  fields  on  integer  time  steps,  n  + 1,  and  update 
the  H  fields  during  half  time  steps,  n  +  Noting  that  the  H  field  cannot  be  updated 
twice  in  one  time  step,  the  last  term  of  Equation  (3.13)  is  approximated  as: 


tt  in 

nx\i,j+l/2,k+l/2 
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After  incorporating  Equation  (3.14)  into  Equation  (3.13)  and  combining  terms: 


tt  in+1/2  _  p.  yt  in— 1/2 

rlx\i,j+l/2,k+l/2  ~  1Ja  '  rLx\i,j+l/2,k+l/2 
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Ez 
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(3.15) 


where  the  coefficients,  Da  and  11;,,  are  composed  of  the  material  parameters  of  the  cell 
being  updated: 


2  ai,j,k ^ 

2^ i,j,k  + 


(3.16) 
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(3.17) 


Analyzing  Equation  (3.15)  in  conjunction  with  Figure  3.1,  it  is  clear  that  the  value  of 
the  Hx  field  component  at  a  particular  cell  is  based  on  the  previous  value  of  that  cell’s  Hx 
field  as  well  as  the  adjacent  orthogonal  E  fields  that  were  updated  one  half  of  a  time  step 
ago.  A  similar  corresponding  analysis  is  true  of  the  other  five  field  components.  Equations 
for  the  other  components  follow  in  much  the  same  manner: 
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where: 
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(3.23) 
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5.^  Absorbing  Boundary  Condition  (ABC) 

Simulations  seek  to  explore  a  model  in  known  conditions.  For  this  simulation,  as 
is  commonly  the  case,  the  antenna  was  to  be  simulated  in  unbounded  free  space.  Since 
it  is  clearly  not  possible  to  simulate  a  model  in  an  environment  whose  extent  is  infinite 
in  all  directions,  the  antenna  was  enclosed  in  an  Absorbing  Boundary  Condition  (ABC). 
An  ABC  seeks  to  absorb  outward  propagating  waves  before  they  can  be  reflected  inward 
toward  the  model. 

In  the  early  days  of  FDTD,  ABC’s  consisted  of  an  expansion  of  the  wave  equation, 
called  a  radiation  operator  [5].  Many  radiation  operators  were  developed  in  the  1970’s 
and  1980’s.  Commonly  used,  was  the  Bayliss-Turkel  operator,  in  which  a  weighted  sum 
of  the  spatial  derivative  in  the  outgoing  direction,  the  spatial  derivative  in  the  transverse 
direction,  and  the  time  derivative  were  taken  of  neighboring  fields  [39].  This  method 
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worked  based  on  a  diminishing  remainder  term.  In  1981,  Mur  proposed  a  finite-difference 
ABC  based  on  the  one-way  approximation  of  the  wave  equation  [5,27].  This  method  was 
used  extensively  for  the  next  15  years  despite  the  drawback  that  Mur’s  absorbtion  was  very 
sensitive  to  frequency.  Schemes  for  surrounding  the  computational  domain  with  a  lossy 
medium  whose  impedance  matches  that  of  free  space  were  attempted  by  several,  however, 
all  had  the  result  of  no  reflection  only  for  plane  waves  at  normal  incidence  [5].  Plane  waves 
at  normal  incidence  can  necessitate  a  rather  large  computational  domain. 

In  1994,  Berenger  made  Mur’s  ABC  all  but  obsolete  when  he  introduced  an  al¬ 
ternative  called  Perfectly  Matched  Layer  (PML)  [33].  PML  is  composed  of  a  layer  cells 
modelling  a  dissipative  material  surrounding  the  FDTD  computational  domain  whose  wave 
impedance  is  perfectly  matched  to  the  space  it  surrounds.  Berenger’s  hypothetical  material 
was  based  on  a  mathematical  model  of  an  anisotropic  medium  in  which  each  component 
of  magnetic  field  is  split  into  two  new  components.  For  example,  in  the  2-D  TE  case, 
Equation  (3.25): 


Mo 


dHz  TT  dEx  dEv 
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dt 
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is  replaced  by  Equations  (3.26)  and  (3.27). 
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v±±zy  .  tt  _  (o 

Mo  Qt  +  vmyHzy  —  (3.27) 

Berenger  demonstrated  that  the  effectiveness  of  this  split-field  PML  was  independent  of 
frequency  and  independent  of  incidence  angles,  providing  the  incidence  angle  is  within  75° 
from  normal  [5].  These  results  for  a  2-D  implementation  were  verified  by  [1].  Katz  also 
verified  Berenger  and  was  the  first  to  extend  PML  to  3-D  [16].  Kantartzis  compiled  a 
comprehensive  comparison  of  PML,  Mur,  higher  order  Mur,  and  higher-order  radiation 
operators,  and  concluded  that  PAIL  was  absolutely  superior  to  all  other  ABC’s  [14],  Also 
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Figure  3.2:  UPML  covering  a  PEC  wall  normal  to  x. 

in  1994,  Chew  and  Weedon  introduced  a  PML  based  on  complex  coordinate  stretching 
variables  whose  effect  is  much  like  Berenger’s  [8,33]. 

Sacks  et  al  developed  a  diagonally  anisotropic  layer,  later  termed  Uniaxial  PML 
(UPML)  by  Taflove,  which  results  in  equivalent  performance  to  split-field  and  stretched 
coordinate  PML’s.  UPML  is  easier  to  implement  since  it  does  not  involve  modifying 
Maxwell’s  equations  and  is  able  to  match  anisotropic  media.  UPML  seeks  to  attenuate 
the  fields  by  creating  an  imaginary  component  of  the  permittivity  in  the  direction  normal 
to  a  PEC  boundary  wall.  Figure  3.2  illustrates  UPML  on  a  wall  that  is  normal  to  the  x 
direction.  Adding  an  imaginary  component  to  ex  will  attenuate  Ex,  and  in  turn,  the  other 
five  field  components.  Permittivity  in  the  PML  in  Figure  3.2  is  chosen  to  be: 

%,-  =  %■(!- <3-28) 

The  reflection  coefficient  for  a  wave  entering  this  UPML  of  a  thickness,  d ,  at  an  arbitrary 
incident  angle,  9 ,  is  given  by: 
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(3.29) 


r(0) 
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Reflection  is  least  when  the  wave  impacts  at  normal  incidence  and  approaches  total 
reflection  as  incidence  approaches  grazing  angle.  Of  course,  UPML  is  intended  to  have 
zero  reflection.  This  appears  to  be  possible  by  simply  making  crx  very  large,  however,  as 
ax  becomes  large,  the  UPML  ceases  to  be  matched  to  the  material  it  is  surrounding  and 
creates  reflection  at  the  surface.  Most  PML’s,  therefore,  use  a  grading  scheme  to  make 
ax  small  at  the  material  boundary  gradually  increasing  to  a  larger  value  at  the  PEC.  The 
UPML  used  in  this  study  is  polynomial  graded,  using: 


(3.30) 


where  ax  is  zero  at  the  material  boundary  and  rises  as  an  mth  degree  polynomial  to  a 
maximum  value  of  (Jx  max  at  the  PEC  wall.  (Tx.max  is  not  without  bound.  If  the  rise  of  loss 
is  too  rapid,  discretization  error  instabilities  will  dominate  any  effects  of  loss.  Equation  3.29 
can  be  updated  with  Equation  3.30  to  become: 
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(3.31) 


Taflove  [39]  states  that  much  empirical  evidence  indicates  that  the  optimum  balance  be¬ 
tween  absorbtion  rate  and  stability  occurs  when  ax,max  is  given  by: 


0.8(m  +  1) 

7/o  dx 


(3.32) 
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Using  Equation  (3.32)  with  Equation  (3.31)  results  in  nearly  140  dB  reduction  at  normal 
incidence  for  a  ten-cell-thick  PML: 


r(o)  =  e^1  (°'K1))  wdcos(0) 

=  e"L6s  =  e~16  =  -139  dB  (3.33) 

The  UPML  derivations  for  the  faces  normal  to  the  other  five  directions  are  analogous. 

3. 3  Materials  specification 

Each  cell  is  a  homogeneous  material  specified  by  it’s  constitutive  parameters:  permit¬ 
tivity,  e.  permeability,  /a,  electrical  conductance,  a,  and  magnetic  conductance,  a*.  These 
parameters  can  be  complex. 

The  conductor  can  be  specified  using  cells  with  the  actual  constitutive  parameters 
of  the  material,  however,  the  conductor  can  instead  be  modelled  as  an  infinitesimally 
thin  layer  of  PEC  by  simply  setting  the  tangential  electric  fields  to  zero  at  the  desired 
cell  boundaries.  The  boundary  condition  on  the  normal  magnetic  field,  h  ■  H  =  0,  is 
automatically  satisfied  by  the  finite  difference  equations  if  the  tangential  electric  fields 
have  been  zeroed,  h  x  E  =  0. 

Figure  3.3  illustrates  a  simulated  PEC  strip  2  cells  in  the  y  direction  by  3  cells  in  the 
x  direction.  Sheen  reported  that  this  method  can  be  used  to  quite  accurately  represent  a 
number  of  simple  microstrip  circuits  and  antennas  [35]. 

In  his  dissertation,  Sheen  later  discussed  a  slightly  more  accurate  method  of  modelling 
PEC,  shown  in  Figure  3.4  [36].  The  difference  equations  developed  in  this  chapter  relied  on 
a  central  difference  approximation,  in  which  case,  PEC  is  equivalent  to  forcing  the  average 
value  of  h  x  E  to  zero. 
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Figure  3.3:  PEC  can  be  modelled  by  setting  the  appropriate  boundary  Electric  field 

components  to  zero  [35]. 
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Figure  3.4:  (a)  PEC  borders  can  be  defined  as  the  last  tangential  E  components  that 

are  zeroed,  (b)  A  more  accurate  method  defines  the  PEC  borders  \  cell  beyond  the  zeroed 
components. 
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3-4  Stability 

The  difference  equations  of  (3.15),  (3.18),  (3.19),  (3. 20), (3. 21),  and  (3.22)  show  that, 
in  this  formulation,  the  fields  are  dependent  on  only  the  fields  in  the  directly  neighboring 
cells.  For  a  stable  solution,  a  propagating  wave,  therefore,  must  not  travel  more  than  one 
cell  in  any  direction  during  one  time  step.  The  maximum  time  step  that  will  prevent  this 
instability  is  given  by  the  dispersion  relation  for  a  plane  wave  propagating  in  a  FDTD  grid. 
To  derive  the  finite  difference  dispersion  relation,  first  reduce  the  source-free  Maxwell’s 
Equations  to  the  wave  equations  for  free  space: 


V2E  = 


1  <92E 

c2  dt2 


V2H  = 


1  <92H 

c2  dt 2 


which  can  be  represented  as: 


(3.34) 

(3.35) 


d 2  d'2  92  \  1  d2il> 

dx2  +  dy2  +  d^)'llJm^~dW  (3'36) 

where  is  any  vector  component  of  E  or  H.  In  general,  the  plane  wave  solution  to  Equa¬ 
tion  (3.36)  is  of  the  form: 


iP(x,  y,  z,  t )  =  Re  ^eo{^kxx-kyy-kzz)  J  (3.37) 

where  kx,  ky,  and  kz  are  the  wavenumber  components  in  rectangular  coordinates.  Equa¬ 
tion  (3.37)  can  be  discretized  in  space  and  time  as: 


5,  nt)  =  Re  {< 


tj  (uJTlt kx  Vr;  x  k y  '^  y  y  z  z  ) 


(3.38) 


where  nx,  ny,  and  nz  are  the  number  of  cells  in  each  of  the  rectangular  directions  and  nt  is 
the  number  of  time  steps,  A t-  Each  cell  is  uniform  and  has  the  dimensions  Ax  x  Ay  x  Az. 
Defining  the  central  difference  operator  in  the  x  direction,  5X,  as: 
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4i/>  = 


i)(nx  +  l,ny,nz,nt)  -  ip(nx  -  \,ny,nz,nt) 


A, 


(3.39) 


and  using  Equation  (3.38)  with  Equation  (3.39)  results  in: 


Sxip  = 


_A  L.  A X  A  U  ^ X 

g  J  K'X  2  _  gJ  2 

aZ 


4 


a7s,"(—  ^ 


(3.40) 


The  second  derivative  operator  is  simply  the  square  of  the  first  derivative  operator: 


=  (44)  '  (4Y0  =  sin2  (3.41) 

Defining  5y,  e>2,  and  <52  analogously,  the  scalar  wave  equation  of  Equation  (3.36)  becomes: 


(&x  + 


sin' 
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2  kyAy 
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1  sin2 
c2  A2 


(3.42) 


which  is  the  dispersion  relation  for  a  plane  wave  propagating  on  the  FDTD  computational 
grid.  Solving  for  lo,  which  must  be  purely  real  for  a  stable  solution,  Equation  (3.42) 
becomes: 
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sin"1^) 


(3.43) 


where 


37 


(3.44) 


£  =  cAt 
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The  upper  bound  of  £  gives  rise  to  the  Courant  condition  for  the  3D  case  of  uniform  cell 
size  free  space  [39]: 


A,  <  \  ,  (3.45) 

cvA  +  ij  +  i? 

If  no  free  space  is  used  in  the  grid,  c  should  be  replaced  in  Equation  (3.45)  with  vrnax , 
which  is  the  maximum  propagation  velocity  through  any  medium  in  the  grid  space. 

The  Courant  condition  of  Equation  3.45,  which  is  also  known  as  the  Courant- 
Fredericks-Levy  (CFL)  criteria,  only  holds  for  the  case  of  this  central-difference  formu¬ 
lation.  Higher  order  difference  equations  may  depend  on  fields  of  more  than  one  cell’s 
distance,  which  allows  a  wave  to  travel  more  than  one  cell  per  time  step  in  a  stable  solu¬ 
tion.  Fewer  time  steps  is  a  great  benefit  of  higher  order  FDTD  formulations  [19]. 

Instability  can  present  itself  in  several  ways  depending  on  the  output  product  of  the 
simulation.  Figure  3.5  shows  the  effects  of  an  instability  that  develops  while  a  traveling 
wave  distribution  is  developing.  Due  to  an  unstable  condition,  the  waveform  develops  small 
spikes,  which  are  barely  visible  near  cell  number  600  of  Figure  3.5(a).  The  instability  may 
or  may  not  originate  in  the  vicinity  of  cell  600  since  the  spikes  could  have  travelled  many 
cells  before  they  were  visible.  The  magnitude  of  the  spikes  increases  as  the  spikes  travel  to 
the  right  in  Figure  3.5(b).  Figure  3.5(c)  shows  the  spikes  impact  the  PML  near  cell  1200 
and  are  partially  attenuated.  The  reduced  magnitude  spikes  reflect  and  begin  traveling  to 
the  left,  as  seen  in  Figure  3.5(d).  The  spikes  again  increase  in  magnitude  (Figure  3.5(e)) 
until  they  overtake  the  waveform  (Figure  3.5(f)).  A  video  of  the  fields  shows  checkerboard 
pattern  of  neighboring  highs  and  lows  spreading  outward  in  all  directions  from  the  point 
of  instability,  as  seen  in  the  stills  of  Figure  3.6.  Notice  that  the  checkerboard  pattern 
alternates  each  timestep. 
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Figure  3.5:  (a)  An  instability  causes  barely  visible  spikes  near  cell  number  600.  The 

magnitude  of  the  spikes  increases  as  the  spikes  travel  to  the  right  in  (b).  (c)  shows  the 
spikes  impact  the  PML  near  cell  1200  and  are  partially  attenuated.  The  reduced  magnitude 
spikes  reflect  and  begin  traveling  to  the  left,  as  seen  in  (d).  The  spikes  again  increase  in 
magnitude  (e)  until  they  overtake  the  waveform  (f). 
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Figure  3.6:  An  instability  quickly  overtakes  the  entire  computational  space. 
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3. 5  Dispersion 

The  first  step  of  applying  an  FDTD  formulation  to  a  structure  is  gridding  the  com¬ 
putational  space.  The  cell  size  chosen  has  implications  for  accuracy  of  the  resulting  data, 
time  required,  and  computing  resources  needed.  As  with  nearly  all  areas  of  CEM,  larger 
discretizations  bring  less  computations,  but  more  error. 

The  most  significant  error  in  the  case  of  FDTD  is  grid  dispersion,  which  is  a  cumu¬ 
lative  phase  error.  Dispersion  is  a  product  of  the  discretization  of  space  as  well  as  time. 
Since  a  wave  travels  a  certain  distance  specified  by  Equation  (3.45)  each  time  step,  a  phase 
front  traveling  along  the  direction  diagonal  to  the  grid  will  propagate  faster  than  the  same 
wave  traveling  along  a  grid  axis.  Simplifying  Equation  (3.42)  for  a  2-D  square  grid  with  a 
monochromatic  traveling  wave: 


cA, 


sm 


A  Sm 


-sm 


2 


(3.46) 


where  A  is  the  length  of  a  grid  cell  edge  and  At  is  the  period  of  a  time  step.  To  satisfy 
the  Courant  condition  of  Equation  (3.45),  let: 


A 

cAt 


=  2 


(3.47) 


Using  the  resolution,  N,  to  equal  the  number  of  cells  per  wavelength,  Equation  (3.46)  can 
be  rewritten  as: 


2  /  7T  \  ,2  /A/ccos0^  .  2  /  Ak  sin  0 


4si,r  Kin)  = sm  v  2 


+  sin 


(3.48) 


Equation  (3.48)  is  a  transcendental  equation  that  can  be  solved  for  k  using  Newton’s 
method  with  the  help  of  trigonometric  identities  [39]: 
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Figure  3.7:  The  magnitude  of  phase  velocity  error,  or  dispersion  [39]. 
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where  m  is  simply  the  iteration  number.  The  normalized  phase  velocity  can  be  found  using 
km  from  the  final  iteration  of  Equation  (3.49)  as: 


vp  2ir 
c  km 


(3.50) 


Figure  3.7  is  a  plot  of  —  showing  that  the  amount  phase  velocity  error  is  proportional  to 
the  resolution  and  the  direction  of  travel.  As  expected,  the  phase  velocity  is  fastest  in  the 
diagonal  direction,  45°,  and  is  slowest  along  a  grid  axis.  This  anisotropy  of  phase  velocity 
is  reduced  as  the  number  of  cells  per  wavelength  is  increased.  Similar  to  MoM,  A/10  is 
a  benchmark  cell  size  that  works  acceptably  in  many  situations.  A/20  or  finer  resolution 
may  be  needed  if  the  scattering  from  a  complicated  structure  is  desired.  Some  researchers 
have  found  cells  as  large  as  A/4  acceptable,  particularly  if  the  waves  modelled  have  few 
interactions  [39]. 
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Figure  3.8:  Rectangular  gridding  results  in  stair-step  approximations  of  curves  and 

diagonal  lines. 

3.6  Geometric  Distortion 

It  doesn’t  take  long  to  notice  the  staircase  effect  from  the  use  of  rectangular  cells  to 
simulate  a  curve,  as  seen  in  Figure  3.8.  The  obvious  solution  is  the  use  of  finer  gridding, 
but  this  leads  to  longer  run  times  and  more  computational  resources.  Curves  are  better 
represented  with  cylindrical  or  spherical  coordinates,  however  the  difference  equations, 
particularly  involving  the  PML,  become  cumbersome.  Conformal  grids  have  been  used 
to  reduce  distortion  of  complex  shapes,  particularly  microstrip  with  curvature  [15,38]. 
Subgridding,  or  finer  gridding  in  only  certain  regions,  has  been  used  by  many  researchers 
beginning  with  Yee  [38].  Attention  must  be  paid  to  the  boundaries  between  unlike  cell 
sizes  to  properly  model  the  field  coupling  between  fine  and  course  grids.  Updating  the 
fields  on  the  boundary  can  be  accomplished  using  interpolation  of  both  space  and  time. 
Kunz  developed  a  method  that  requires  two  runs,  the  first  to  update  the  fields  in  the 
course  grid,  the  second  to  update  the  fields  in  the  fine  grid  using  the  coarse-grid  fields  as  a 
boundary  [19].  A  common  solution  to  geometric  distortion  is  a  hybrid  technique  employing 
FDTD  and  MoM  or  FDTD  and  Geometrical  Theory  of  Diffraction  (GTD). 

3. 7  Source 

There  are  two  types  of  sources  that  can  be  used:  a  hard  source  or  a  soft  source.  For 
a  hard  source,  the  field  at  one  or  more  cells  is  forced  to  have  a  certain  value  regardless  of 
the  neighboring  fields.  Waves  cannot  travel  through  these  hard  source  cells  and  appear  to 
reflect  off  in  the  same  manner  as  if  the  cell  was  a  perfect  electric  conductor  (PEC).  For  a 
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Figure  3.9:  Cubic  growth  over  the  first  few  periods  eliminates  problematic  high  frequency 
components  of  rapid  transitions. 

soft  source,  the  introduced  value  is  only  added  to  the  existing  field  of  the  cell  or  cells.  Soft 
sources  are  beneficial  when  you  do  not  have  prior  knowledge  of  the  fields  in  the  vicinity 
of  the  source  as  these  cells  allow  waves  to  travel  through.  Hard  sources  have  the  ability  of 
being  turned  off  more  easily  than  soft  sources. 

The  waveform  of  the  source  also  must  be  considered.  A  sharp  start  or  stop  transition 
in  the  source  wave  will  create  messy  propagating  waves.  These  discontinuities  contain  many 
higher  frequency  components  that  can  disturb  the  surrounding  fields  to  such  an  extent  that 
a  source  may  even  enter  runaway  mode  and  fail  to  turn  off.  Taflove  notes  that  FDTD  has 
an  intrinsic  low  pass  filter  effect  since  computations  are  only  based  on  neighboring  cells 
[39] .  Theoretically,  the  high  frequency  components,  if  left  long  enough,  would  eventually 
disappear.  Typically,  these  discontinuities  are  eliminated  by  a  gaussian  pulse,  which  is  a 
sinusoidal  function  attenuated  by  an  exponential,  as  seen  in  Equation  (3.51). 


sourcent 


sm(untAt)  ■  e*- 


(3.51) 


where  £  is  an  attenuation  constant  that  is  typically  related  to  the  resolution,  or  number 
of  cells  per  wavelength.  For  this  study,  I  chose  a  constant  source,  so  I  was  not  concerned 
with  a  stop  transition.  I  decided  to  use  a  sinusoidal  source  with  a  cubic  growth  over  the 
first  three  periods,  which  is  illustrated  in  Figure  3.9. 
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IV.  Simulation  Development 


This  chapter  presents  the  methodology  by  which  an  FDTD  simulation  was  developed 
to  analyze  antennas  of  different  geometries. 

4.1  Hagness-Willis  Code 

To  save  time,  I  obtained  a  working  3-D  FDTD  program  with  a  UPML  absorbing  layer. 
The  UPML  was  a  10-cell-thick,  4th  order  polynomial  graded,  and  PEC-backed  on  all  six 
faces.  This  “beta”  version  Matlab  code  was  provided  by  Dr.  Susan  Hagness  and  one  of  her 
Ph.D.  students,  Keely  Willis,  both  of  the  University  of  Wisconsin  Computational  Electro¬ 
magnetic  Laboratory.  The  Hagness-Willis  code  will  be  featured  in  the  upcoming  edition  of 
Taflove  and  Hagness’s  book  Computational  Electromagnetics:  The  Finite- Difference  Time 
Domain  Method,  3rd  Ed.  The  program  was  extensively  modified  for  my  purposes  and,  as 
a  result,  I  used  little  more  than  their  UPML  implementation. 

The  first  design  simulated  was  the  220  mm  long  full  width  Thiele  antenna  using  a 
solid  wall  down  the  centerline,  as  seen  in  Figures  4.1  and  4.2.  Since  the  UPML  of  the 
Hagness-Willis  code  was  matched  to  free  space,  I  surrounded  the  entire  antenna  with  a 
20-cell-thick  layer  of  free  space.  Above  the  antenna,  I  included  an  additional  20  cells,  for  a 
total  of  40  cells  of  free  space  between  the  top  UPML  boundary  and  the  top  of  the  antenna. 
The  open  substrate  was  extended  outward  one  antenna-width  in  both  y  directions.  The 
following  paragraphs  detail  the  procedure  I  used  to  pare  away  unnecessary  cells  revealing 
only  the  minimum  number  of  cells  required  to  accurately  model  the  structure. 

4-2  Copper 

I  began  modelling  the  copper  as  1-cell-thick  layer  of  £o>  Mo>  cr  =  5.8  x  107  S/m. 
Designs  were  fabricated  with  Rogers  Corporation  5870  RT/Duroid  with  1  ounce  copper 
cladding,  which  is  35/mr  thick.  If  the  FDTD  cells  were  nearly  cubic,  the  substrate  would 
need  to  be  22  cells  thick  and  the  entire  computational  space  would  require  nearly  900  mil¬ 
lion  cells.  From  smaller  trials  that  I  have  run,  I  estimate  that  this  would  call  for  upwards  of 
1  TB  of  memory  for  even  the  lossless  case.  I  abandoned  this  method  and  instead  approxi- 
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Figure  4.1:  The  y  —  z  cross-section  slice  of  the  Thiele  Full  Width  (TFW)  antenna 

surrounded  by  free  space  (not  to  scale).  The  source  cell  is  shown  in  gold.  The  x  —  z  slice 
of  Figure  4.2  is  shaded  in  blue. 
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Figure  4.2:  The  x  —  z  longitudinal-section  slice  of  the  TFW  antenna  surrounded  by  free 
space  (not  to  scale).  The  source  cell  is  shown  in  gold.  The  y  —  z  slice  of  Figure  4.1  is 
shaded  in  blue. 
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mated  the  copper  conductor  as  a  perfect  electric  conductor  (PEC)  of  zero  thickness.  This 
approximation  allowed  for  much  larger  cells,  decreasing  the  total  number.  For  example,  by 
making  the  substrate  only  5  cells  thick,  the  simulation  is  cut  from  900  million  cells  down 
to  roughly  50  million  cells,  a  94  %  decrease. 

4.3  PML 

The  UPML  of  the  Hagness-Willis  code  is  matched  to  free  space,  therefore,  my  antenna 
was  first  surrounded  by  a  layer  of  free  space,  which  was  then  surrounded  by  the  UPML 
cells.  To  extract  the  propagation  constant  of  the  forward  traveling  wave,  I  needed  to 
eliminate  reflections  from  the  free  space  boundaries  at  the  ends  of  the  substrate.  This 
was  accomplished  by  extending  the  substrate  directly  into  the  UPML,  and  modifying 
the  affected  UPML  cells  to  match  the  substrate.  The  resulting  grid  space  is  seen  in 
Figures  4.3  and  4.4.  Notice  that  the  UPML  layer  under  the  ground  plate  has  also  been 
removed  since  it  is  redundant.  The  UPML  now  absorbs  all  outward  propagating  waves  in 
the  substrate  allowing  the  forward  traveling  wave  to  develop  exclusively.  At  this  point,  the 
UPML  is  inhomogeneous  in  the  z  direction.  If  the  electric  flux  density,  D^,  is  discontinuous 
between  the  two  materials,  Gauss’s  Law  states  that  a  surface  charge  will  develop  at  the 
boundary.  A  charge  building  between  two  bordering  cells  is  a  source  of  instability  for  the 
FDTD  method.  This  problem  was  averted  in  the  UPML  development,  since  the  imaginary 
component  of  the  UPML  permittivity  in  Equation  (3.28)  was  chosen  to  be  normalized 
to  free  space  permittivity  instead  of  the  material’s  permittivity.  The  number  of  cells  to 
simulate  the  antenna  after  the  removal  of  the  surrounding  free  space  was  reduced  over  71% 
to  14.2  million  cells,  which  was  still  much  too  large. 

4-4  Material  Surrounding  the  Antenna 

Since  the  PML  absorbs  nearly  all  of  the  outward  propagating  energy,  the  amount  of 
free  space  above  and  the  amount  of  open  substrate  extending  outward  from  the  structure 
has  no  effect  on  kx  as  long  as  the  PML  Is  at  least  two  cells  away  from  the  structure.  The 
resulting  grid  space  with  only  the  required  number  of  cells,  a  further  reduction  of  almost 
80  %,  is  seen  in  Figures  4.5  and  4.6. 
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Figure  4.3:  The  cross-section  of  the  TFW  antenna  extending  into  the  UPML  (not  to 

scale).  The  source  cell  is  shown  in  gold.  The  x  —  z  slice  of  Figure  4.4  is  shaded  in  blue. 


Figure  4.4:  The  longitudinal-section  of  the  TFW  antenna  extending  into  the  UPML 

(not  to  scale).  The  source  cells  are  shown  in  gold.  The  y  —  z  slice  of  Figure  4.3  is  shaded 
in  blue. 
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Figure  4.5:  The  y  —  z  slice  of  the  TFW  antenna  extending  into  the  UPML  with  unneeded 
material  removed  (not  to  scale).  The  source  cells  are  shown  in  gold.  The  x  —  z  slice  of 
Figure  4.6  is  shaded  in  blue. 
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Figure  4.6:  The  x—z  slice  of  the  TFW  antenna  extending  into  the  UPML  with  unneeded 
material  removed  (not  to  scale).  The  source  cells  are  shown  in  gold.  The  y  —  z  slice  of 
Figure  4.5  is  shaded  in  blue. 
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4-5  Resolution  and  Cell  Size 


I  ran  several  trials  to  determine  if  the  cell  size  could  be  changed.  By  keeping  the 
cross  section  square  and  elongating  the  cells  in  the  x  direction,  Figure  4.7(a)  shows  that  at 
the  leaky  band’s  center  frequency,  7  GHz,  there  is  less  than  1%  error  of  (3X  if  the  cell  size  is 
6:1:1  or  less.  The  error  is  with  respect  to  the  transverse  resonance  solution.  Figure  4.7(b) 
is  a  plot  of  the  error  due  to  varying  cell  sizes  while  keeping  the  cells  square  in  the  x  and  y 
directions.  Cell  size  must  be  1.5:1. 5:1  or  less  to  stay  within  0.5%  of  transverse  resonance. 
Figure  4.7(c)  shows  that  keeping  the  height  of  the  substrate  to  at  least  5  cells  ensures 
agreement  with  transverse  resonance  to  within  0.5%.  Using  a  cell  size  of  5:1:1  with  a 
five-cell-thick  substrate,  the  grid  of  Figures  4.5  and  4.6  is  only  572,000  cells. 

The  worst  resolution  will  occur  at  areas  with  the  shortest  wavelength  and  largest  cells. 
All  subsequent  trials  used  five-cell-thick  substrate  and  sizes  of  either  3:1:1  or  5:1:1.  Since 
the  only  substrate  thickness  used  was  787  yrn,  the  cells  were  472.2  x  157.4  x  157.4  yrn 
or  787  x  157.4  x  157.4  ym,  respectively.  The  shortest  wavelength  encountered  in  the 
simulation  was  in  the  substrate  outside  of  the  antenna,  given  by: 


substrate 


(4.1) 


Therefore,  the  worst  resolution  is  given  by  Figure  4.8  for  the  largest  dimension  of  each 
of  the  two  cell  sizes.  The  resolution  of  the  traveling  wave  inside  the  structure  was  in  the 
hundreds  of  cells  per  wavelength. 


4-6  Source 

Many  source  geometries  and  locations  were  tried  to  find  the  optimum.  At  least 
one  source  cell  must  be  more  than  two  cells  from  the  PML  or  too  much  energy  will  be 
absorbed  for  the  traveling  wave  to  adequately  develop.  A  source  cell  centered  vertically 
in  the  substrate  allows  a  traveling  wave  of  the  largest  magnitude  to  develop.  A  source 
cell  directly  on  the  antenna’s  open  edge  prevents  establishment  of  the  traveling  wave  by 
allowing  the  source  energy  to  propagate  away  from  the  antenna  and  not  in  the  structure. 
Following  these  findings,  all  contiguous  geometries  consisting  of  single  cells,  rows  of  cells 
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3.5 


(a) 


(b) 


(c) 


Figure  4.7:  (a)  Plot  of  the  error  of  (3X  as  a  function  of  elongating  the  x  direction  of  the 

cells  while  keeping  the  other  dimensions  square,  (b)  Plot  of  the  error  of  /3X  as  a  function  of 
cell  size  for  cells  square  in  the  x  and  y  directions,  (c)  Plot  of  the  error  of  Bx  as  a  function 
of  the  height  of  the  substrate  in  cells.  For  all  figures,  the  error  is  with  respect  to  the 
transverse  resonance  solution  and  the  curve  is  a  cubic  least-squares  fit  of  the  data  points. 
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(a)  (b) 

Figure  4.8:  The  worst  resolution  encountered  in  Thiele  antenna  simulations  with  5-cell- 
thick  substrate  for  cell  sizes  of  (a)  3:1:1  and  (b)  5:1:1. 

in  each  of  the  three  rectangular  directions,  and  walls  of  cells  normal  to  each  of  the  three 
directions,  performed  identically.  All  subsequent  trials  used  a  single  source  cell  midway 
between  the  ground  plate  and  the  conductor  strip,  two  cells  inside  from  the  open  edge  of 
the  antenna  and  five  cells  from  the  PML. 

4-1  Loss 

Simulating  substrate  with  loss  requires  nearly  twice  the  computer  RAM  compared 
to  the  lossless  case  since  each  field  component  of  each  cell  must  have  a  real  and  imaginary 
portion.  Rogers  Duroid  5870  high  frequency  laminate  has  a  loss  tangent  of  only  0.0012. 
Neglecting  this  loss  showed  no  noticeable  effect  on  the  extracted  propagation  constant,  as 
can  be  seen  in  Figure  4.9.  All  subsequent  trials  were  treated  as  lossless  freeing  50%  of  the 
memory  allowing  a  larger  simulation. 

4-8  Precision 

Matlab  defaults  to  double  precision  unless  specified  otherwise.  Trials  run  in  single 
precision  used  roughly  40%  less  memory  and  the  results  were  consistent  with  double  pre¬ 
cision,  as  seen  in  Figure  4.10.  All  subsequent  trials  remained  double  precision  unless  the 
size  of  the  simulation  required  more  memory  than  was  available. 
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Figure  4.9:  No  noticeable  error  of  the  extracted  propagation  constant  due  to  neglecting 
loss  for  the  Thiele  Half  Width  (THW)  antenna. 


Figure  4.10:  Propagation  constant  extracted  from  the  fields  following  a  single  precision 
simulation  does  not  differ  significantly  from  the  double  precision  results. 
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Figure  4.11:  The  propagation  constant  data  extracted  from  the  THW  antenna  is  indis¬ 
tinguishable  from  the  TFW  antenna. 

4-9  Full  Width  vs.  Half  Width 

Figure  4.11  illustrates  the  agreement  between  the  propagation  constants  of  the  Thiele 
Half  Width  (THW)  and  Thiele  Full  Width  antennas  (TFW)  to  within  1%  over  the  entire 
leaky  band.  Simulation  of  the  THW  antenna  uses  approximately  45%  less  memory  than 
the  TFW  simulation. 

4-10  Determination  of  Leakage  Constant  and  Phase  Constant 

The  objective  of  the  antenna  simulations  was  limited  to  providing  the  propagation 
constant  of  the  vertical  component  of  the  electric  field,  E2,  inside  the  substrate  between  the 
top  conductor  and  the  bottom  ground  plate.  The  E2  data  was  retrieved  from  a  single  row  of 
cells  stretching  the  length  of  the  antenna,  as  depicted  in  green  in  Figure  4.12.  As  discussed 
in  Chapter  2,  the  attenuation  constant,  ax,  and  the  phase  constant,  /3X,  components  of 
the  propagation  constant  give  an  accurate  means  of  comparing  the  bandwidth,  leakage 
rate,  main  beam  direction,  and  approximate  far-field  pattern  of  different  traveling  wave 
antennas. 
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Figure  4.12:  4.12(a)  is  a  vertical  slice  and  4.12(b)  is  a  horizontal  slice  of  the  compu¬ 

tational  domain  (not  to  scale).  Ez  was  retrieved  from  the  cells  in  green  to  be  used  to 
determine  the  propagation  constant.  The  source  cell  is  shown  in  gold. 
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Figure  4.13:  The  raw  Ez  data  retrieved  following  simulation. 

Typical  raw  data  following  a  simulation  is  plotted  in  Figure  4.13.  The  source  is  seen  as 
the  large  spike  at  cell  number  13.  At  cell  420,  the  PML  begins  to  dampen  the  waveform.  At 
least  two  periods  of  the  waveform  of  Figure  4.13,  starting  from  the  first  absolute  maximum 
(cell  72),  were  normalized  for  analysis.  A  recursive,  least-square  procedure  was  used  to 
best  fit  a  known  exponential  curve.  The  peak  values  were  compared  to  find  ax  and  the  zero 
crossing  locations  were  compared  to  find  (3X-  Figure  4.14  shows  the  best  fit  curve  on  top  of 
the  FDTD  data.  The  procedure  was  automated  to  process  an  entire  set  of  simulation  trials 
and  output  the  propagation  constants  to  an  excel  spreadsheet.  This  procedure  works  well 
for  most  frequencies,  however,  its  use  for  the  lowest  15%  of  a  leaky  band  is  problematic 
for  two  reasons. 

First,  good  agreement  with  a  best  fit  method  is  not  adequate  with  less  than  two 
periods  of  data.  The  wavelength  of  the  traveling  wave  at  the  lower  end  of  the  leaky  band 
is  more  than  five  times  the  wavelength  at  the  highest  frequency  in  the  band,  as  seen  in 
Figure  4.15.  Therefore,  the  computational  domain  must  be  five  times  longer,  which  requires 
five  times  the  memory.  Even  worse,  the  processing  time  will  be  increased  by  more  than 
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X-direction  cell  number 


Figure  4.14:  The  data  from  Figure  4.13  has  been  normalized  and  matched  with  its 

best-fit  exponential  curve. 


Figure  4.15:  The  phase  constant  wavelength,  Xp,  is  much  longer  for  lower  frequencies. 
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25  times,  since  each  time  step  will  take  five  times  longer  and  the  waveform  will  need  five 
times  the  steps  to  propagate  the  length  of  the  antenna.  A  large  simulation  is  difficult  and 
time-consuming,  but  possible. 

The  second  problem  was  a  bigger  roadblock.  For  frequencies  at  which  ^  is  greater 
than  0.1,  the  waveform  attenuates  too  rapidly  to  find  sufficient  peaks  or  zero  crossings. 
Instead,  a  much  simpler  method  was  developed  from  the  logarithm  of  the  normalized 
“FDTD”  waveform  of  Figure  4.14  using  Equation  (4.2): 


Ez  oc  e(a*~j0*)x 

In  Ez  oc  lne^*-^*^  =  axx  —  j(3xx  (4.2) 

As  shown  in  Figure  4.16,  ax  is  simply  the  slope  of  the  peaks  of  the  In  Ez  waveform,  shown 
as  a  dotted  line,  and  j3x  is  found  from  the  separation  of  nulls  using: 

0X  =  ^  (4.3) 

A/3 

If  the  simulation  has  run  long  enough  to  ensure  that  the  traveling  wave  distribution  has  a 
constant  wavelength,  only  ^  is  needed  for  determination  of  the  propagation  constant.  A 
computational  space  that  is  just  one-half  wavelength  long  is  a  |  reduction  from  the  best-fit 
method  that  required  three  wavelengths. 

4-11  Validating  FDTD  Code 

4-11.1  Transverse  Resonance.  The  FDTD  simulation  was  validated  by  com¬ 
parison  of  the  extracted  ax  and  (3X  with  the  propagation  constant  given  by  a  transverse 
resonance  solution,  which  was  verified  by  a  comparison  with  the  Steepest  Descent  Contour 
(SDC)  method  by  Lee  and  Oliner  [20,30].  Figure  4.17  shows  a  transmission  line  model 
that  is  applicable  to  the  cross  section  of  the  Menzel,  TFW,  and  THW  antennas.  Each 
structure  can  be  modelled  as  a  dielectric-filled  parallel  plate  waveguide  of  admittance,  Yq£, 
terminated  at  one  end  by  a  short  circuit  and  the  other  end  by  admittance  T).  The  E  null 
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Figure  4.16:  The  natural  logarithm  of  the  simulation  data  can  be  used  to  determine  the 
propagation  constant. 

generated  in  the  EHi  mode  by  a  vertical  wall  or  transverse  slots  is  represented  by  a  short 
circuit.  Yt  is  an  approximation  of  the  admittance  of  an  open  edge  of  microstrip  developed 
by  Chang  and  Kuester  [7, 17]  using  the  Weiner-Hopf  technique  to  analyze  a  TEM  wave 
that  is  completely  reflected. 

Chang  and  Kuester  state  applicability  to  only  thin  substrates  in  which: 


h  <C  — (4.4) 

W  y/£jl 

For  a  substrate  thickness,  h,  of  787 nm.  Equation  (4.4)  requires  frequencies  <C  40  GHz, 
which  is  five  times  higher  than  the  highest  frequency  of  this  leaky  band.  At  6.7  GHz, 
Figure  4.18  shows  that  FDTD  and  transverse  resonance  begin  to  disagree  as  the  height 
of  the  substrate  increases  past  the  thin  criteria,  near  h=  1.1  mm.  By  Equation  4.4,  the 
frequency  should  be  much  less  than  approximately  28.4  GHz  at  this  point.  This  suggests 
that,  for  the  THW  antenna  (er  =  2.33  and  w=  7.5  mm),  Equation  (4.4)  is  more  precisely: 


h  < 


28.4  4.2 

2vr6.7  y/ejlo  lu  ^/ejlo 


(4.5) 
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l/2  width 


]  Short 
Circuit 


Transverse 


(a)  (b) 

Figure  4.17:  (a)  is  a  transmission  line  circuit  that  approximates  the  cross  section  of  each 

of  the  three  structures  in  (b)  operating  in  mode  EHi.  The  cross  sections  in  (b)  represent 
Menzel  (top),  TFW  (center),  and  THW  (bottom). 


Figure  4.18:  The  effect  of  the  height  of  the  substrate  on  the  propagation  constant  at 

6.7  GHz  for  the  THW  antenna  (er  =  2.33;  w  =  7.5  mm.)  The  dashed  line  is  a  quadratic 
least  square  of  the  FDTD  data  points. 
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The  transverse  resonance  relation,  Equation  (4.6),  must  hold  for  all  points  in  the 
transverse  direction,  y. 


r right (y)  '  D left{y )  —  1 


(4.6) 


The  reflection  coefficient  from  the  admittance  of  the  end  of  the  microstrip,  Y),  is 
unity  with  a  phase  shift,  x  [7, 17].  At  a  point  y  =  ya  just  to  the  right  of  Yt. 


r  rightly  a)  =  ~e~j2k^ 

(4.7) 

I"  leftiVa)  =  eJX 

(4.8) 

where  k  is  the  wavenumber  in  the  substrate  and  w  is  the  width  of  the  structure.  Equa¬ 
tion  (4.6)  becomes: 


_e-f2fcf  .  eix 
X  —  kw 


1 

±n7r  n  =  1, 3,  5, ... 


(4.9) 


For  the  EHi  mode,  n=l. 

The  approximation  can  be  modified  to  meet  changes  in  h,  w,  and  £>,  providing  that 
Equation  (4.4)  is  still  valid.  It  cannot,  however,  be  applied  to  curves  or  tapers  since 
the  Yt  approximation  is  no  longer  valid.  Figure  4.19  shows  that  the  transverse  resonance 
approximation  is  in  agreement  within  1%  of  the  FDTD-derived  (3X  data  over  the  entire 
leaky  band. 

4-11.2  Measurements.  Validation  was  also  shown  by  comparison  with  existing 
measurements.  The  Thiele  Half  Width  antenna  was  fabricated  at  the  Radiation  and  Scat¬ 
tering  Compact  Antenna  Laboratory  (RASCAL)  by  Dr.  Thiele  using  vias  spaced  A/ 10  to 
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Figure  4.19:  The  FDTD  simulation  of  the  THW  is  in  agreement  with  the  transverse 

resonance  approximation. 

form  the  wall.  Both  far-held  pattern  measurements  as  well  as  near-held  probing  measure¬ 
ments  were  conducted. 

The  near  held  measurements  produced  the  data  seen  in  Figure  4.20,  from  which  the 
propagation  constant  was  extracted  in  much  the  same  manner  as  with  the  FDTD  data. 
Figure  4.21  shows  that  the  data  is  not  very  smooth,  but  /3X  corresponding  to  the  peaks 
does  agree  to  within  3%  of  both  FDTD  and  transverse  resonance. 

Far-held  H-plane  measurements  were  taken  at  RASCAL’s  compact  range.  The  set¬ 
up  is  depicted  in  Figure  4.22.  Figure  4.23  shows  the  normalized  far-held  H-plane  E  held 
pattern  of  the  THW  antenna  made  with  vias  by  Dr.  Thiele.  For  comparison,  the  data  is 
plotted  along  with  the  line  source  pattern  of  the  corresponding  data  from  both  the  FDTD 
simulation  and  the  transverse  resonance  approximation.  The  magnitude  of  the  backward 
lobe  is  higher  for  the  measured  THW  antenna  demonstrating  a  lower  attenuation  constant, 
ax.  The  location  of  the  measured  main  lobe  is  closer  to  endhre,  which  indicates  a  higher 
(3X  than  the  numerical  solutions. 
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Figure  4.20:  The  raw  data  resulting  from  probing  the  near  field  of  the  Thiele  half  width 
antenna  at  6.7  GHz. 


Near.mat 


Figure  4.21:  Extraction  of  the  propagation  constant  from  the  near  field  measurements 

(6.7  GHz). 
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(b) 

Figure  4.22:  Far-field  H-plane  measurements  of  the  E  field  taken  in  RASCAL’s  compact 
range.  (a)Horizontal  view  and  (b)  View  from  above  showing  placement  of  antenna  under 
test  and  the  probe  of  the  receive  horn  (left). 


The  estimated  propagation  constant  of  the  Thiele  antenna  can  be  determined  by  a 
best  fit  with  the  line  source.  Figure  4.24  shows  the  line  source  pattern  matched  to  that  of 
the  THW  antenna.  As  discussed  in  Chapter  2,  the  main  beam  position  is  only  within  1°  if 
the  antenna  is  over  5  long.  The  length  of  all  antennas  fabricated  was  limited  to  19  cm 
due  to  material  availability.  This  length  is  less  than  3  \p  at  6.7  GHz  and  less  than  4 
at  7.2  GHz.  To  extract  the  propagation  constant,  5%  was  added  to  the  line  source  ax  and 
/3X  to  accurately  match  the  main  beam  location  within  1%.  The  propagation  constants  of 
the  THW  antenna  determined  from  Figure  4.24  are  summarized  in  Table  4.1. 


Table  4.1: 


Line  source  best  fit  of  the  Thie 


.e  antenna  made  with  vias. 


6.7  GHz 

7.2  GHz 

Measured  kxj ko 

0.715  -  j0.019 

0.85  -  j0.0135 

FDTD  kx/ko 

0.677  -  j0.040 

0.816  -  j0.0252 

(3X  Relative  error 

5.7% 

4.2  % 
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6.7  GHz 


(a) 


7.2  GHz 


Angle  from  endfire  (degrees) 

(b) 

Figure  4.23:  Far-field  radiation  pattern  of  the  THW  fabricated  with  vias  compared  to 
the  line  source  pattern  generated  with  FDTD  and  Transverse  Resonance  data  (a)  at  6.7 
GHz  and  (b)  at  7.2  GHz. 


65 


6.7  GHz 


7.2  GHz 


(a) 


(b) 


Figure  4.24:  Measured  data  with  a  best  fit  line  source  to  estimate  the  propagation 

constant  (a)  at  6.7  GHz  and  (b)  at  7.2  GHz. 


4-12  Curvature 

The  effect  of  curvature  on  the  propagation  constant  was  tested  by  constructing  a 
180°  bend  of  the  THW  design.  Figure  4.25  illustrates  the  view  from  above  the  curved 
computational  domain,  looking  in  the  +z  direction.  The  substrate  is  shown  as  cyan,  the 
PML  is  maroon,  and  the  top  conducting  strip  is  colored  green.  For  viewing,  the  wall  has 
been  colored  blue  and  the  cells  from  which  were  recorded  are  colored  orange.  The 
source  is  the  single  dark  cell  near  (20,140).  The  antenna  was  bent  with  the  wall  on  the 
inside  (Figure  4.25(a))  as  well  as  with  the  wall  on  the  outside  (Figure  4.25(b)). 

Geometrical  distortion  was  expected  to  pose  a  problem.  Figure  4.26  is  a  blown-up 
view  of  the  cells  at  a  curve.  A  cell  size  of  1.5:1. 5:1  was  used  to  keep  the  cells  square  in  the 
plane  of  curvature  to  minimize  distortion. 

The  most  difficult  task  involved  in  simulating  curvature  in  a  rectangular  coordinate 
system  is  accurately  processing  the  data.  As  seen  in  Figure  4.26,  the  cells  from  which  the 
field  data  is  taken  are  not  on  a  precise  curve.  Worse,  the  cell  data  is  not  stored  in  the 
order  the  traveling  wave  propagates.  These  matters  were  solved  by  storing  the  angle  to 
each  held  cell  with  respect  to  the  lower  PML  boundary  as  seen  by  the  cell  at  the  center 
of  the  semi-circle,  cell  (138,141).  The  output  data  from  this  simulation  was,  therefore,  a 
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(a) 


(b) 

Figure  4.25:  The  computational  space  of  the  THW  antenna  curved  180°  in  Matlab.  (a) 
shows  the  wall  (in  blue)  on  the  inside  while  (b)  shows  the  wall  on  the  outside 
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Figure  4.26:  The  curved  THW  antenna  is  approximated  by  rectangular  cells. 

vector  called  EZDATA,  which  contained  the  component  field  data  from  the  orange  cells 
in  Figure  4.4.25,  and  a  vector  of  EZDATA’ s  corresponding  angles,  named  ANGLES. 

Figure  4.27  depicts  the  processing  using  the  angle  information.  Figure  4.27(a)  shows 
the  ANGLES  data  following  simulation.  The  phase  is  unwrapped  in  Figure  4.27(b).  Next, 
both  EZDATA  and  ANGLES  are  re-ordered  based  on  the  ANGLES  data  to  put  the  data 
in  the  order  encountered  by  the  traveling  wave  propagating  in  the  antenna,  as  seen  in 
Figure  4.27(c).  Finally,  the  vertical  axis  of  Figure  4.27(d)  shows  the  ANGLES  have  been 
translated  into  arc  distances  from  the  source  cell  to  allow  computation  of  the  propagation 
constant  from  the  EZDATA.  The  propagation  constant  can  then  be  extracted  in  the  same 
manner  as  the  straight  antennas. 
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(a)  (b) 


(c)  (d) 


Figure  4.27:  Following  the  curved  THW  simulation,  processing  of  the  (a)  raw  angle  data 
involved  (b)  unwrapping  the  phase,  (c)  ordering,  and  (d)  converting  the  angle  to  an  arc 
distance  from  the  source. 
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V.  Results 


Analysis  of  the  Thiele  Half  Width  (THW)  antenna  began  by  a  examining  the  field  dis¬ 
tribution  of  its  predecessors,  the  Menzel  antenna  and  the  Thiele  Full  Width  (TFW) 
antenna.  One  of  the  primary  uses  of  determining  the  propagation  constant  of  a  traveling 
wave  antenna  is  bandwidth  prediction.  Antennas  were  simulated  with  varying  geometrical 
parameters  and  the  resulting  effects  on  the  propagation  constants  were  analyzed  to  find 
means  of  increasing  the  bandwidth.  A  preliminary  look  into  the  effects  of  an  array  of 
THW  elements  determined  the  required  spacing  to  prevent  mutual  interaction  of  multiple 
elements.  The  current  method  of  fabricating  the  THW  antenna  is  an  imprecise,  tedious 
process  of  drilling  holes  and  soldering  vias.  An  improved  fabrication  technique  replacing 
the  vias  with  conducting  copper  tape  is  documented. 


5. 1  Reduction  of  Memory  for  Simulation 

The  simulation  development  steps  taken  in  the  previous  section  allowed  the  trials  to 
be  accurately  run  on  personal  computers  with  only  1  GB  of  RAM.  Table  5.1  summarizes 
means  that  can  be  used  to  eliminate  unnecessary  cells  and  the  corresponding  decrease  in 
memory  use  observed.  These  steps  reduced  the  initial  simulation  of  900  million  cells  over 
99.9%  to  only  572,000.  Unless  otherwise  annotated,  all  trials  modelled  the  THW  antenna 
as  a  five-cell-thick  lossless  substrate  with  uniform  cells  of  size  3:1:1  using  double  precision. 
These  approximations  enable  simulations  of  up  to  1  million  cells. 


Table  5.1:  Ways  the  size  of  the  simulation  was  reduced. 


Assumption 

Memory  Reduction 

Modelling  copper  as  zero  thickness 

over  94  % 

Extraction  of  kx  using  log 

over  80  % 

Cell  size  5:1:1  (3:1:1) 

80  %  (60  %) 

Elimination  of  unneeded  substrate 

nearly  80  % 

5-cell-thick  substrate 

nearly  80  % 

No  free  space  surrounding  structure 

over  70  % 

Lossless  substrate 

50  % 

Remove  non-excited  side  of  TFW 

approx  45  % 

Single  precision 

approx  40  % 

Total 

over  99.9  % 
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5.2  Menzel  antenna 


The  wall  of  the  Thiele  design  is  clearly  superior  to  transverse  slots  at  suppression  of 
the  fundamental  mode.  The  transverse  slots  of  the  Menzel  design  created  field  distributions 
that  made  FDTD  extraction  of  the  propagation  constant  difficult  for  higher  frequencies 
and  virtually  impossible  for  the  lower  half  of  the  leaky  band,  as  shown  in  Figure  5.1.  The 
patterns  of  Figures  5.1(b)  and  5.1(c)  suggest  destructive  interference  from  the  presence 
of  the  EHo  mode  in  addition  to  EHi.  The  right  half  of  the  pattern  enclosed  in  a  dotted 
box  of  Figures  5.1(d)  and  5.1(e)  is  a  waveform  useable  for  propagation  constant  extraction 
and  yields  data  within  5%  of  both  the  transverse  resonance  approximation  and  the  THW 
simulation.  An  attempt  to  improve  the  waveform  by  enlarging  the  transverse  slots  is  seen 
in  Figure  5.2.  The  lower  frequencies  are  still  unusable,  but  the  usable  portion  of  the  higher 
frequencies  has  been  increased  slightly. 

5.3  Thiele  Full  Width  (TFW)  antenna 

Figure  5.3  illustrates  an  investigation  into  the  fields  inside  the  TFW  antenna.  The 
source  cell  seen  in  the  upper  left  corner  of  Figure  5.3(a)  has  an  Ez  component  with  a 
negative  (red)  value  that  induces  positive  (blue)  fields  in  the  immediately  surrounding 
cells  on  the  excited  (top)  side.  As  the  positive  (blue)  fields  £11  the  excited  side,  negative 
(red)  fields  are  induced  into  the  non-excited  side  (bottom)  past  the  vertical  PEC  wall. 
Figures  5.3(b)-5.3(e)  show  that  these  induced  fields  form  a  traveling  wave  in  the  non- 
excited  side  that  travels  in  the  same  +x  direction  and  has  roughly  the  same  wavelength, 
\p,  as  the  excited  side.  Since  the  energy  propagated  along  a  longer  path,  the  non-excited 
side  is  not  180°  out  of  phase  but  less  than  90°.  Like  those  of  the  excited  side,  the  fields 
of  the  non-excited  side  attenuate  as  they  travel.  This  attenuation  indicates  that  radiation 
is  also  occurring  on  the  non-excited  side.  Due  to  the  less  than  90°  phase  difference  with 
respect  to  the  excited  side,  the  radiation  from  the  non-excited  side  is  likely  decreasing  the 
ability  of  the  structure  to  radiate.  Extracting  kx  from  simulations  does  not  support  any 
difference  between  the  THW  and  the  TFW.  However,  Dr.  Thiele  reports  that  previous 
measurements  comparing  the  far-field  pattern  of  both  antennas  indicate  that  the  THW 
does  in  fact  radiate  with  higher  gain  than  the  TFW. 
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View  from  above  (+z)  at  6.7  GHz 
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Figure  5.1:  (a)  The  original  Menzel  antenna  as  simulated  in  Matlab  at  6.7  GHz.  The 

PML  appears  maroon,  PEC  green,  and  substrate  blue.  The  antenna  yielded  the  following 
traveling  wave  distributions:  (b)  at  6.2  GHz,  (c)  at  6.7  GHz,  (d)  at  7.2  GHz,  and  (e)  at 
7.7  GHz.  The  only  useable  data  is  surrounded  by  a  dotted  line. 
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View  from  above  (+z)  at  6.7  GHz 
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Figure  5.2:  (a)  The  Menzel  antenna  with  larger  transverse  slots  as  simulated  in  Matlab 

at  6.7  GHz.  The  PML  appears  maroon,  PEC  green,  and  substrate  blue.  This  antenna 
yielded  the  following  traveling  wave  distributions:  (b)  at  6.2  GHz,  (c)  at  6.7  GHz,  (d)  at 
7.2  GHz,  and  (e)  at  7.7  GHz.  The  only  useable  data  is  surrounded  by  a  dotted  line. 
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(c) 


Top  view:  Ez,  time  step  =  1920 
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Figure  5.3:  The  excited  side  (top  of  each  figure)  of  the  TFW  antenna  induces  a  traveling 
wave  on  the  non-excited  side  (bottom)  that  also  travels  to  the  right  (+x)  and  attenuates 
as  it  radiates.  The  source  cell  can  be  seen  in  the  upper  left  corner  of  (a),  (c),  and  (d). 
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Figure  5.4:  The  propagation  constant’s  dependence  on  dielectric  constant  at  6.7  GHz 

for  the  THW  antenna  (h  =  0.787  mm;  w  =  7.5  mm.) 

5-4  Modifying  dimensions  to  meet  bandwidth  specifications 

A  bandwidth  centered  around  a  desired  center  frequency,  /c,  can  be  achieved  by 
scaling  the  width  of  the  conducting  strip,  the  height  (or  thickness)  of  the  substrate,  and/or 
the  permittivity  of  the  substrate.  The  bandwidth  can  be  increased,  to  a  limited  extent,  by 
the  selection  of  the  substrate  material.  As  a  fraction  of  the  center  frequency,  the  bandwidth 
will  be  unaffected  by  altering  the  height  and  width,  although  fc  can  be  shifted  readily. 
The  choice  of  substrate  material  and  thickness  is  usually  dictated  by  cost  or  availability  of 
material,  therefore,  the  width  is  the  easiest  of  the  three  to  manipulate. 

5.4.1  Varying  Dielectric  Constant.  Figure  5.4  illustrates  the  relationship  between 
the  relative  permittivity  of  the  substrate  and  the  propagation  constant.  Figure  5.5  shows 
the  bandwidth  as  a  function  of  substrate  permittivity  overlayed  with  common  substrates. 
Bandwidth  increases  rapidly  as  the  substrate  dielectric  constant  nears  that  of  free  space. 
The  drawback  of  lower  dielectric  constant  is  a  very  low  ax  across  most  of  the  leaky  band, 
as  seen  in  Figure  5.6.  Low  ax  results  in  little  energy  radiating  per  unit  length. 
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Figure  5.5:  The  bandwidth  as  a  function  of  substrate  permittivity  for  the  THW  Antenna 
(h  =  0.787  mm;  w  =  7.5  mm.)  Common,  commercially  available  substrates  are  overlayed 


TRh787w150e133 


Figure  5.6:  The  normalized  ax  and  (3X  curves  in  the  leaky  band  for  the  THW  antenna 
using  substrate  with  er=1.33  (h  =  0.787  mm;  w  =  7.5  mm.) 


76 


Figure  5.7:  The  effect  of  the  width  of  the  conductor  on  the  propagation  constant  at 

6.7  GHz  for  the  THW  antenna  [er  =  2.33;  h  =  0.787  mm.) 

5-4-2  Varying  Height.  Figure  4.18  illustrates  the  relationship  between  the  height, 
or  thickness,  of  the  substrate  and  the  propagation  constant.  Like  permittivity,  the  height 
of  the  substrate  is  usually  dictated  by  the  material  available.  All  antennas  simulated  and 
fabricated  for  this  work  used  material  that  was  0.787  mm  thick.  As  mentioned  in  the 
previous  section,  the  transverse  resonance  solution  includes  an  approximation  that  limits 
its  applicability  to  thin  substrates.  The  difference  in  j3x  between  FDTD  and  transverse 
resonance  becomes  noticeable  for  heights  greater  than  approximately  1.1  mm. 

5-4-3  Varying  Width.  Figure  5.7  shows  that  the  propagation  constant  is  very 
sensitive  to  the  width  of  the  conductor  strip.  As  little  as  0.1  mm  difference  in  width 
will  cause  as  much  as  10%  error  in  f3x .  This  sensitivity  to  width  questions  the  ability  to 
fabricate  the  wall  with  vias,  since  placement  to  even  as  little  as  0.1  mm  is  difficult. 

5-4-4  Frequency  Scaling.  The  leaky  frequency  band  can  be  scaled  up  or  down 
by  scaling  the  width  and  height  inversely  while  using  the  same  permittivity.  For  example, 
Table  5.2  shows  the  effect  on  frequency  by  halving  the  height  and  width.  While  the 
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Figure  5.8:  The  center  frequency  has  been  scaled  up  by  a  factor  of  four  by  using  a 

quarter  of  the  height  and  a  quarter  of  the  width. 


bandwidth  appears  to  double,  the  effective  bandwidth  remains  the  same  percentage  of  the 
center  frequency,  fc. 

Table  5.2:  Scaling  the  frequency  by  a  factor  of  two. 


Bandwidth 

2.4  GHz 

4.8  GHz 

fc 

7.08  GHz 

14.16  GHz 

W 

15  mm 

7.5  mm 

h 

787  fj,m 

393.5  firri 

£r 

2.33 

2.33 

Figure  5.8  shows  the  effect  of  using  one- fourth  the  width  and  one- fourth  the  height  of 
the  THW.  The  center  frequency  has  quadrupled  to  28.32  GHz.  Frequency  scaling  does 
not  prove  useful  for  reduction  of  the  FDTD  simulation  since  cross  section  ratio  of  height 
to  width  is  unchanged.  However,  this  property  could  be  useful  to  shrink  the  antenna  to 
meet  measurement  facility  size  constraints.  Conversely,  fabrication  may  be  made  easier  by 
increasing  the  size  of  the  antenna. 
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Figure  5.9:  The  THW  antenna  curved  90°  with  a  radius  of  approximately  93  nun. 

5. 5  Curvature 

Constant  curvature  with  the  open  end  on  the  outward  side  was  simulated  at  radii 
of  3.36,  4.23,  5.32,  and  9.3  cm.  Only  90°  was  possible  to  simulate  for  the  largest  radius, 
while  180°  was  used  for  all  others.  The  size  of  the  computational  domain  necessitated  only 
single  precision.  Not  enough  information  is  available  to  extract  kx  for  frequencies  whose 
\/3  is  greater  than  twice  the  length  of  the  arc.  Figure  4.15  shows  the  relation  between 
frequency  and  the  wavelength  of  /3X.  The  curvature  trials  are  summarized  in  Table  5.3. 
Figure  5.9  shows  that  curvature  increases  bandwidth  by  flattening  (3X,  predominantly  for 
the  lower  frequencies,  while  keeping  ax  relatively  unaffected.  Figures  5.10  through  5.12 
indicate  that  as  curvature  decreases,  the  bandwidth  decreases. 


Table  5.3:  Summary  of  curvature  trials. 


180u 

180u 

180u 

90u 

Radius  (cm) 

3.36 

4.23 

5.32 

9.3 

Arc  length  (cm) 

10.6 

13.3 

16.7 

14.6 

Largest  A g  possible  (cm) 

21.2 

26.6 

33.4 

18.6 

Lowest  /  possible 

6.2  GHz 

6.1  GHz 

6.0  GHz 

6.1  GHz 

Approx  Bandwidth 

2.8  GHz 

2.7  GHz 

2.6  GHz 

2.5  GHz 
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Figure  5.10:  The  THW  antenna  curved  180°  with  a  radius  of  approximately  53 


Figure  5.11:  The  THW  antenna  curved  180°  with  a  radius  of  approximately  42 


mm. 


mm. 
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Figure  5.12:  The  THW  antenna  curved  180°  with  a  radius  of  approximately  34  mm. 

Curvature  with  the  wall  on  the  outside  appears  to  hamper  the  ability  of  the  antenna 
to  set-up  the  EHi  mode.  As  seen  in  Figure  5.13,  which  is  data  from  a  180°  arc  of  radius 
4.23  cm,  destructive  interference  indicates  the  presence  of  the  EHo  mode. 

5.6  Multiple  Elements 

From  Equation  (2.16)  it  is  clear  that  the  main  beam  is  steerable  in  the  longitudinal 
direction  from  near  endfire  to  near  broadside.  A  linear  array  of  these  elements  should, 
therefore,  be  able  to  scan  in  two  dimensions.  A  first  step  to  developing  such  an  array 
is  to  determine  the  required  spacing  between  neighboring  elements.  Two  elements  were 
simulated  next  to  each  other  at  7.2  GHz  using  double  precision.  The  results  were  nearly 
identical  regardless  of  whether  the  second  element  was  excited  or  not.  The  effect  on 
propagation  constant  by  a  neighboring  element  is  shown  in  Figure  5.14.  There  appears  to 
be  little  interaction  between  elements  if  they  are  separated  by  at  least  0.25AS  and  virtually 
no  interaction  if  the  spacing  is  over  0.4AS,  where: 

As  =  -y=  (5.1) 

V  £r 
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(a)  (b) 

Figure  5.13:  Curvature  with  the  wall  on  the  outside  produces  field  distributions,  seen 

here  at  (a)6.4  GHz  and  (b)  7.4  GHz,  that  suggest  the  presence  of  multiple  modes. 


(a)  (b) 

Figure  5.14:  (a)  Actual  error  of  \oix\  and  (b)  relative  error  of  f3x  from  placing  another 

element  a  fraction  of  a  wavelength  from  the  open  end  of  the  antenna.  For  both  figures,  the 
error  is  with  respect  to  the  propagation  constant  of  a  single  element. 
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Figure  5.15:  The  propagation  constant  from  FDTD  simulations  is  unaffected  by  the 

removal  of  substrate  outside  of  the  wall. 

5.7  Simplified  Fabrication 

Fabricating  the  THW  antenna  using  vias  soldered  to  the  conducting  strip  and  ground 
plate  is  time  consuming  and  inaccurate.  Precise  placement  of  a  0.52  mm  thick  wire  is 
quite  difficult.  As  previously  mentioned,  Figure  5.7  shows  that  an  error  in  width  of  as 
little  as  0.1  mm  can  affect  the  propagation  constant  by  as  much  as  10%.  The  spacing 
of  the  vias  may  also  create  a  reactive  component  to  the  wall  which  could  further  alter 
the  propagation  constant.  Replacing  the  vias  with  conductive  copper  tape  was  quick  and 
simple.  Figure  5.15  does  not  support  the  need  for  the  substrate  outside  the  wall,  therefore, 
the  wall  could  be  created  by  cutting  the  substrate  and  ground  plate  away  up  to  the  location 
of  the  wall,  applying  copper  tape  to  form  the  wall,  and  then  soldering  a  replacement  ground 
plate.  The  resulting  antenna  is  pictured  in  Figure  5.16. 

Figures  5.17  and  5.18  are  comparisons  of  far-field  patterns  among  the  two  THW 
antennas  and  the  line  source  pattern  from  the  FDTD-produced  propagation  constant. 
The  magnitude  of  the  backward  lobes  indicates  that  attenuation  constant,  ax,  of  the 
THW  antenna  made  with  conductive  tape  more  closely  matches  the  FDTD  data  than  does 
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Figure  5.16:  THW  antenna  fabricated  without  vias. 


the  antenna  utilizing  vias.  Position  of  the  main  lobe  of  each  antenna  indicates  that  /3X  is 
nearly  identical  for  each. 

Table  5.4  summarizes  the  estimated  error  of  the  propagation  constant  of  the  THW 
antenna  fabricated  with  conducting  copper  tape.  Propagation  constant  was  estimated  with 
a  line  source  pattern  best  fit  match  of  the  main  and  backward  lobes.  At  6.7  GHz,  fJx  has 
slightly  more  error  than  the  antenna  made  with  vias.  Fabrication  error  making  the  width 
of  the  conducting  strip  slightly  wider  than  designed  is  attributed.  The  width  of  the  strip 
is  more  crucial  and  leads  to  more  error  for  lower  frequencies  than  higher  frequencies. 


Table  5.4:  The  THW  antenna  made  with  copper  tape. 


6.7  GHz 

7.2  GHz 

Measured  kx/ko 

0.732  -  j0.034 

0.86  -  j0.023 

FDTD  kx/ko 

0.677  -  j0.040 

0.816  -  j0.0252 

Px/ko  Relative  error 

8.2  % 

5.4  % 
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6.7  GHz 


Figure  5.17:  Comparison  at  6.7  GHz  of  THW  using  vias,  THW  using  conductive  tape, 
and  the  line  source  -  FDTD  pattern  estimate. 


7.2  GHz 


Figure  5.18:  Comparison  at  7.2  GHz  of  THW  using  vias,  THW  using  conductive  tape, 
and  the  line  source  -  FDTD  pattern  estimate. 
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VI.  Conclusions 


This  thesis  provides  the  groundwork  that  will  enable  development  of  a  lightweight, 
inexpensive,  aerodynamic,  and  broadband  antenna.  Whether  for  radar  or  communi¬ 
cation,  an  antenna  with  these  properties  would  be  a  force  multiplier  for  the  smaller,  limited 
payload  air  vehicles  the  United  States  Air  Force  will  pursue  in  the  coming  years. 

6. 1  Summary  of  Results 

A  three  dimensional  Finite  Difference  Time  Domain  (FDTD)  simulation  was  con¬ 
structed  to  simulate  various  leaky  wave  microstrip  antennas.  FDTD  was  shown  to  be  able 
to  extract  the  propagation  constant,  which  is  useful  in  predicting  antenna  performance 
measures  such  as  bandwidth  and  main  bean  direction.  The  antenna  geometry  was  easily 
manipulated  in  FDTD  to  meet  various  testing  requirements  of  different  structures.  Sev¬ 
eral  means  were  found  to  decrease  the  number  of  cells  required  to  accurately  simulate  the 
antennas  allowing  elaborate  simulations  with  only  a  standard  desktop  PC  with  1  GB  of 
memory.  The  availability  of  high  performance  computing,  which  is  typically  required  for 
such  modelling,  was  shown  to  be  not  necessary  for  most  geometries. 

The  FDTD  simulation  provided  results  within  1%  of  a  transverse  resonance  approxi¬ 
mation  of  the  same  structure.  The  transverse  resonance  model  has  been  validated  by  other 
researchers  with  the  Steepest  Descent  Contour  [20,30].  Measurements  were  also  used  for 
validation.  A  best  fit  line  source  far  field  pattern  was  matched  to  the  measured  pattern 
to  estimate  the  propagation  constant  for  comparison  with  the  simulations.  The  ability  to 
quantitively  compare  simulated  antenna  propagation  constants  with  propagation  constants 
of  measured  antennas  was  previously  not  available. 

The  Thiele  Half  Width  (THW)  antenna  was  shown  to  be  an  improvement  over  both 
the  Thiele  Full  Width  (TFW)  antenna  as  well  as  the  Menzel  antenna.  The  propagation 
constant  of  the  THW  antenna  behaved  equivalently  to  the  TFW  antenna  as  anticipated 
by  image  theory.  In  addition,  the  non-excited  side  of  the  TFW  antenna  was  shown  to 
degrade  the  phase  difference  across  the  structure’s  width  preventing  effective  radiation. 
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The  vertical  wall  of  the  THW  antenna  was  found  to  be  more  capable  of  blocking  the 
fundamental  mode  than  the  antenna  developed  by  Menzel  [23]. 

Two  THW  antenna  elements  were  excited  in  close  proximity  to  test  mutual  inter¬ 
action.  No  effect  to  the  propagation  constant  of  elements  was  observed  if  the  spacing 
distance  was  greater  than  0.4A.  This  minimum  spacing  will  allow  multiple  THW  elements 
to  adequately  radiate  as  an  array. 

Bandwidth  improvements  were  found  to  be  inversely  related  to  curvature  of  the  THW 
antenna.  Geometrical  distortion  may  be  a  large  source  of  error  for  the  curved  antennas 
simulated.  Therefore,  these  findings  need  confirmation  with  measured  or  analytical  results. 
The  far-held  line  source  pattern  used  is  only  valid  for  a  straight  antenna.  A  curved  antenna 
can  be  fabricated  and  measured,  however,  the  propagation  constant  cannot  be  estimated 
with  the  line  source  developed.  A  curved  line  source  would  be  needed  for  comparison  of 
the  propagation  constant  with  simulated  data. 

The  bandwidth  of  the  THW  antenna  can  also  be  increased  by  using  a  lower  permit¬ 
tivity  dielectric  constant;  however,  a  commercially  available  substrate  could  not  be  found 
lower  than  2.2.  Lowering  the  dielectric  constant  has  the  drawback  of  a  very  low  ax  across 
much  of  the  leaky  band,  which  hinders  an  effective  radiator.  Varying  the  width  of  the 
conducting  strip  or  the  height  of  the  substrate  were  found  to  move  the  center  frequency, 
but  did  not  change  the  bandwidth  as  a  percentage  of  center  frequency.  The  width  and 
height  can  be  varied  together  to  scale  the  center  frequency.  Frequency  scaling  can  be  used 
to  decrease  the  length  of  the  antenna  to  meet  measurement  facility  size  constraints  or  to 
increase  the  size  of  the  antenna  for  ease  of  fabrication. 

6.2  Follow-on  Work 

Knowledge  gained  from  this  thesis  about  the  nature  of  the  propagation  constant 
as  a  function  of  curvature  and  element  spacing  is  critical  for  future  efforts  to  produce  a 
spiral  design.  The  Thiele  antenna  shows  promise,  but  more  work  needs  to  be  done  before 
an  effective,  broadband,  broad  pattern  antenna  can  be  made.  Further  research  of  the 
THW  would  be  aided  by  a  faster  simulation  using  more  memory  and  a  transformation  of 
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the  simulated  fields  to  the  far-field.  Areas  of  analysis  that  could  next  be  tackled  include 
continuation  of  the  curved  geometry  investigations  and  exploration  of  a  full  THW  array. 
An  important  topic  not  addressed  in  this  work  is  the  feed. 

6.2.1  FDTD  Simulation.  A  larger  simulation  is  advantageous  for  a  host  of 
reasons.  Of  course,  more  cells  mean  higher  resolution  and  less  geometrical  distortion.  A 
limitation  of  the  FDTD  method  was  that  the  propagation  constant  was  not  able  to  be 
determined  when  the  length  of  the  antenna  was  less  than  half  of  A^.  Simulations  were 
limited  to  900,000  cells  for  double  precision  and  1,600,000  cells  for  single  precision  for 
1  GB  of  computer  memory.  More  memory  would  not  prohibit  extraction  of  the  long  Xg 
propagation  constants  of  low  frequencies.  Instead  of  terminating  the  structure  into  PML, 
a  larger  simulation  would  allow  the  antenna  to  instead  be  surrounded  by  free  space,  which 
would  model  the  fields  and  currents  of  the  actual  antenna  in  addition  to  propagation 
constant  extraction.  The  loss  tangent  of  certain  substrates  may  be  large  enough  to  require 
a  lossy  model,  necessitating  much  more  memory. 

Computer  memory  is  the  factor  limiting  the  simulation  size.  The  solution  is  to  create 
a  parallel  version  of  the  code  to  allow  several  computers  working  the  simulation  simulta¬ 
neously.  The  E  field  components  of  each  cell  are  determined  from  that  cell’s  previous  time 
step  as  well  as  the  directly  adjacent  H  fields.  The  H  field  components  update  analogously. 
The  E  fields  are  not  updated  at  the  same  time  as  the  H  fields,  and  vice  versa.  Therefore, 
breaking  up  the  computational  space  allowing  multiple  processing  of  separate  portions  of 
the  computational  space  simultaneously  is  possible.  Many  methods  of  parallelizing  Matlab 
have  been  developed,  including  MatlabMPI,  which  is  available  in  Wright-Patterson  AFB’s 
High  Performance  Computing  facility. 

The  length  of  time  needed  to  run  the  simulation  was  not  a  problem  due  to  the 
large  number  of  PC’s  that  were  available  for  this  research.  If  time  is  an  issue,  as  would 
likely  be  the  case  if  the  simulations  were  larger,  Matlab  should  be  abandoned  for  the  field 
update  loop  of  the  current  FDTD  program.  Fortran  would  provide  a  vast  improvement  in 
computational  speed. 


A  curved  line  source  approximation  is  needed  to  evaluate  the  curved  simulation 
against  a  pattern  measurement  of  a  fabricated  curved  design.  However,  the  line  source 
approximation  gives  an  indication  of  only  the  main  lobes  of  the  far-held  pattern.  Since 
the  fields  immediately  above  the  simulated  antenna  are  known,  a  near-held  to  far-held 
transformation  would  provide  a  more  exact  simulated  pattern  for  comparison. 

6.2.2  THW  antenna.  The  best  bandwidth  improvements  involve  increasing  ax 
while  not  lowering  /3X.  This  was  shown  to  be  the  effect  of  curvature.  Since  the  bandwidth  of 
the  THW  antenna  is  increased  by  curvature,  the  natural  next  step  is  a  spiral.  Taper  would 
likely  further  increase  bandwidth  and  this  could  be  investigated  with  the  FDTD  simulation. 
A  substrate  with  a  graded  dielectric  constant,  although  quite  difficult  to  construct,  could 
provide  the  benefits  of  increased  bandwidth  and  a  decreased  backward  lobe. 

Since  the  main  beam  is  frequency  steerable,  a  linear  array  of  THW  antennas  could 
allow  a  two  dimensional  scanning  capability.  Limited  investigation  shows  that  mutual 
coupling  does  not  appear  to  be  a  factor  in  the  leaky  region  if  elements  are  separated  by 
0.4A.  Simulation  of  an  array  would  require  much  more  computation  than  even  the  curved 
trials.  Parallel  processing  and  a  near-field  to  far-field  transformation  would  enhance  this 
investigation. 

6.2.3  Feed.  A  proper  impedance  match  with  the  source  is  vital  to  efficient  power 
use  by  an  antenna.  Gain  of  the  antenna  would  also  be  enhanced.  The  simulated  fields  inside 
the  THW  antenna  should  provide  the  information  needed  to  characterize  the  impedance  of 
the  structure.  An  optimum  feed  could  then  be  realized  using  the  required  matching  stubs, 
although  this  match  may  only  be  effective  for  a  narrow  frequency  band.  Since  the  goal  is  a 
broadband  antenna,  a  balun  might  be  necessary  to  provide  the  desired  insertion  loss  over 
the  entire  leaky  band. 

6.2.4  S-Parameter  Measurements.  Sheen  developed  an  S-parameter  extraction 
technique  to  de-embed  alpha  and  beta  from  network  analyzer  measurements  [37].  This 
method  could  be  useful  for  geometries  that  are  difficult  to  simulate.  The  drawback  to  this 
method  is  the  high  cost  of  building  multiple  designs. 
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Appendix  A.  Matlab  Code 


TH.m 


function  TH(f , steps , sing) 


'/,  3-D  FDTD  code  with  UPML  absorbing  boundary  conditions 

'/,  Date  of  this  version:  21  Jan  2005 

•/.  PATTERNED  AFTER  CODE  WRITTEN  BY: 

'/,  Keely  J.  Willis,  Ph.D.  Student 

'/,  UW  Computational  Electromagnetics  Laboratory 

'/,  Director:  Prof.  Susan  C.  Hagness 

/ 

•/.  INPUTS:  ‘f’  IS  THE  EXCITATION  FREQUENCY  IN  (Hz) 

•/.  ‘steps’  BEGINS  TIME  STEPPING  IF  ‘1’,  SKIPS  TIME  STEPPING  TO  ONLY  DISPLAY 

•/.  SIMULATION  PARAMETERS  IF  ‘O’ 

•/.  ‘sing’  SELECTS  SINGLE  PRECISION  IF  ‘1’,  DOUBLE  PRECISION  IF  ‘O’ 

7. 

7. 

'/,  SIMULATION  ‘name’  of  the  form  ‘THexample67 ’  where: 

7. 


•/.  ‘TH’  IS  A  THIELE  HALF  WIDTH  (THW)  ANTENNA 

•/.  or  ‘TF’  IS  A  THIELE  FULL  WIDTH  (TFW)  ANTENNA 

•/.  or  ‘MENZEL’  IS  THE  MENZEL  ANTENNA 

•/.  or  ‘  THCURVEOUT  ’  IS  A  CURVED  THW  WITH  THE  WALL  OUTSIDE 

•/.  or  ‘  THCURVEIN  ’  IS  A  CURVED  THW  WITH  THE  WALL  INSIDE 

•/.  ‘example’  IS  THE  SPECIFIC  CONFIGURATION,  SUCH  AS: 

•/.  ‘loss’  WHICH  INCLUDES  LOSS 

•/.  ‘sing’  WHICH  IS  SINGULAR  PRECISION 

•/.  ‘norm’  WHICH  IS  LOSSLESS  DOUBLE  PRECISION  3:1:1 

•/.  ‘smaller’  WHICH  IS  A  THINNER  ‘norm’ 

•/.  ‘ exxx ’  WHICH  IS  ‘norm’  WITH  ANOTHER  DIELECTRIC 

•/.  ‘ hxxx ’  WHICH  IS  ‘norm’  WITH  ANOTHER  SUBSTRATE  THICKNESS 

•/.  ‘rxx’  WHICH  IS  THE  RADIUS  OF  CURVATURE  (cm) 

•/.  ‘67’  IS  THE  EXCITATION  FREQUENCY  OF  6.7  GHz 
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•/. 

•/. 

•/.  EXAMPLE  OF  THE  SIMULATION  MATLAB  DISPLAY: 


1  - * 

'/,  THexample67  began  at  13-Feb-2005  12:04:53. 

'/,  Computational  grid  is  158790  cells. 

'/,  Resolution  is:  62.24  cells/wavelength. 

'/,  Structure  size:  66.9  X  7.50  X  0.787  (mm) 

'/,  Cell  size  ratio:  2.99  :  1.00  :  1.00 

'/,  Time  steps:  13020  (dt=  0.3437  pS) 

7. 

'/,  FDTD3D  preprocessing  took  1.55  seconds. 

•/. 

'/,  Period  1  of  30  took  9.14  minutes. 

•/.  (13- Jan-2005  12:14:07) 

•/. 

•/. 

•/.  <  DISPLAY  FOR  EACH  PERIOD  OF  SOURCE  > 


y  ii 

y  ii 

'/,  Period  30  of  30  took  9.18  minutes. 

•/.  (13- Jan-2005  17:01:02) 

7. 

•/.  End  at  13-Feb-2005  17:01:35. 

•/. 

•/.  OUTPUTS : 

•/.  ‘THexample67  .mat  ’  :  A  FILE  CONTAINING  ALL  OF  THE  SIMULATION  DATA  THAT  WILL  BE 

•/.  USED  TO  EXTRACT  THE  PROPAGATION  CONSTANT. 

7. 

•/.  ‘THexample67QUART.mat’  :  A  FILES  SAVED  AFTER  A  FOURTH  OF  THE  SIMULATION  IS 

•/.  COMPLETE.  VERY  USEFUL  TO  ESTIMATE  THE  REMAINING  SIMULATION  TIME. 

•/. 

•/.  ‘  MTHexample67 .  mat  ’  :  AN  OPTIONAL  FILE  CONTAINING  FRAMES  OF  THE  FIELDS. 

•/.  (WARNING:  THIS  FILE  CAN  GET  HUGE  "200  MB) 

7. 
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7.  ‘THexample67.avi’  :  AN  OPTIONAL  FILE  OF  MOVIE  IN  COMPRESSED  FORMAT 

7 


tic 


Begin=  datestr  (now)  ;  7.  STORE  BEGINING  TIME  OF  SIMULATION 

f  pr intf  ( ’  \n\n - *  ’ )  ; 

name=  [mf ilename , ’ example ’ ,num2str (f /le8)] ; 

7.  CREATE  A  SIMULATION  NAME 


fprintf  ( ’\n"/„s  began  at  ’/0s  .  ’  ,name ,  Begin)  ; 


curve =  0; 
mov=  0; 
av=  0 ; 
if  av 


7.  ’1’  TO  UPDATE  PEC  FOR  CURVE 
7.  TO  SAVE  A  MOVIE  OF  THE  FIELDS 
7.  ’1’  TO  SAVE  .avi  MOVIE  (ONLY  IF  mov=  ’1’) 
7.  CREATE  A  .avi  FILE  TO  BE  WRITTEN  TO 


aviobj  =  avif ile ( [name avi ’], ’fps 12 quality 100) ; 


end; 


7.  FREQUENCY 

if  sing  c=  single  (299792458)  ;  "/,  speed  of  light  (m/s) 

else  c=  299792458;  end 

omega=  2*pi*f ;  7,  angular  frequency  (rad/s) 

lambda=  c/f;  "/,  wavelength  (m) 

k=  2*pi/lambda;  7>  wavenumber  (1/m) 


7.  MATERIAL  PARAMETERS  :  m=l :  free  space  m=2:  Substrate 


u0=  pi*4e-7 ; 
e0=  c~-2/u0; 
lossTan=  0.00; 

7.  lossTan=  0.0012; 


7.  permeability  of  free  space  (H/m) 
7.  permittivity  of  free  space  (F/m) 
7.  lossless  case 
7.  for  RT/Duroid  5870 
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er=[l  2 . 33+j *lossTan*2 . 33] ;  "/,  relative  permitivity 
e=  eO*er; 

sig=  [zeros  (1 ,  length(er)  )  ]  ;"/0  electrical  conductance 
ur=  [ones (1 , length (er) )  ];  %  relative  permeability 

u=  ur*uO; 

for  m=  1 :  length(er)  '/,  CREATE  COEFFICIENTS  SPECIFYING  MATERIAL  FOR  EACH  CELL 

C . cl (m)=  (2*e (m)-sig(m) *dt)/ (2*e(m)+sig(m) *dt) ; 

C . c2(m)=  2*e (m) *dt/ (2*e (m)+sig(m) *dt)  ; 

C . c3(m)=  C . cl (m) ; 

C.c4(m)=  l/(2*e (m) . ~2+e (m) *sig(m) *dt)  ; 

C.c5(m)=  (2*e(m)+sig(m)*dt) ; 

C. c6(m)=  (2*e(m)-sig(m)*dt) ; 

D. dl(m)=  1; 

D.d2(m)=  dt ; 

D.d3(m)=  1; 

D.d4(m)=  l/(2*e (m) *u0) ; 

D . d5(m)=  2*e (m) ; 

D . d6(m)=  2*e (m) ; 

end 

•/.  GRID  PARAMETERS 

load  Compare 

pABx=round (pABx* 100) / 100 ; 

B=  pB(f ind(pABx==f /le9) ) ; 

L=  c/B/f;  */.  wavelength  predicted  by  Transverse  Resonance 

X=  2*L;  */.  length  (m)  (for  approx  2  wavelengths) 

Y=  (7 . 42e-3)  ;  '/  width  of  strip  (m)  GUESS:  HALF  CELL  SHORTER 

Zs=  787e-6;  '/  thickness  of  substrate  (m) 

dz=  Zs/5;  */.  5  cell  thick  substrate 

if  curve  ’/  1.5: 1.5:1  cell  ratio 
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dx=  X/round(X/(l .  5*dz)  )  ;7o  size  of  <x>  cell  edge  (m) 
dy=  Y/round(Y/(l .  5*dz) )  ;°l  size  of  <y>  cell  edge  (m) 

r=  2/pi*X;  7.  RADIUS  OF  CURVATURE 

jstruct=  round ( (r+Y-r/sqrt (2) ) /dy)  +4; 
istruct=  round(sqrt(2)*r/dx) ; 
if  mod(istruct ,2) 

istruct=  round(sqrt (2) *r/dx) +1 ; 

end 

kstruct=  Zs/dz; 


else  */.  3:1:1  cell  ratio 

istruct=  round(X/3/dz) ;  i  size  of  structure  (#cells)  in  <x> 
jstruct=  round(Y/l/dz)+4;'/,  size  of  structure  (#cells)  in  <y> 
kstruct=  Zs/dz;  7,  size  of  structure  (#cells)  in  <z> 

dx=  X/istruct;  7,  size  of  <x>  cell  edge  (m) 

dy=  .0075/(jstruct-3.5)  ;70  size  of  <y>  cell  edge  (m) 


7.  HALF-CELL  CORRECTION  MADE  TO  WIDTH  OF  CONDUCTING  STRIP 


end 


if  f/le9  >  6.9 
elseif  f/le9  >  6.6 
else 

buffer=  2; 
pb=  pml+buffer; 


pml=  6;  7,  thickness  of  PML  region  (#  cells) 
pml=  8; 

pml=  12;  end; 

7.  free  space  buffer  above  structure 


it=  istruct+2*pml ;  7> 
jt=  j struct+2*pml ;  7, 
kt=  kstruct+pb;  7o 
numcells=  (it*jt*kt;  7o 


f pr intf ( ’ \nComputational  grid 


size  of  total  sim  (#cells) 
size  of  total  sim  (#cells) 
size  of  total  sim  (#cells) 
TOTAL  NUMBER  OF  CELLS 
is  7. .  Of  cells  .’ ,numcells)  ; 


in  <x> 
in  <x> 
in  <x> 


7.  CELL  SIZE  /  TIME  STEP 

DX=  1/dx;  DY=  1/dy;  DZ=  l/dz;7.  MULTIPLICATION  IS  SLIGHTLY  FASTER  THAN  DIVISION 
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res=  lambda/ sqrt (2 . 33) /max( [dx  dy  dz] ) ;  7.  THE  COARSEST  RESOLUTION  IN  THE  GRID 

fprintf  ( ’  \nResolution  is:  "/, .  2f  cells/wavelength.’,  res); 
fprintf  ( ’\nStructure  size:  "/,.lf  X  '/, .  2f  X  "/0.3f  (mm)’,... 

istruct*dx*1000,  (jstruct-3.5)*dy*1000,  kstruct*dz*1000) ; 

cellx=  dx/min([dx  dy  dz] ) ; 
celly=  dy/min([dx  dy  dz] ) ; 
cellz=  dz/min([dx  dy  dz] ) ; 

fprintf  (’ \nCell  size  ratio:  "/,. 2f  :  "/, .  2f  :  '/, .  2f  cellx,  celly  ,  cellz  ); 

dt=  . 9/ (c*sqrt (DX~2  +  DY~2  +  DZ~2)); 

7,  ENSURE  dt  IS  LESS  THAN  COURANT  CRITERIA 
periods=  30;  "/,  #  periods  of  source  (more  for  lower  freqs) 

TimeSteps=  round(periods*round(l/(f *dt) ) ) ;  7,  TOTAL  NUMBER  OF  TIME  STEPS 

fprintf  ( ’\nTime  steps:  "/,.0f  (dt=  '/,.4f  pS) TimeSteps, dt*lel2) ; 


•/.  SOURCE 

s0=  100;  7  SOURCE  MAGNITUDE 

n=  3;  7  number  of  periods  of  ramp 

numramp=  n*round(l/ (f *dt) ) ; 

ramp=  ones (1 ,TimeSteps) ; 

ramp(l  :numramp)=  ( (1  :numramp) /numramp)  .  ~3 ;  "/,  CUBIC  RAMP  OVER  FIRST  THREE  PERIODS 

ramp=  ramp(l : TimeSteps) ; 

S0=  sO*ramp; 

source=  SO  .  *sin(omega*  (  (1 :  TimeSteps)  )  *dt)  ;  7.  CONSTANT  FREQUENCY  SOURCE 

clear  ramp  sO  SO; 


7.  INITIALIZE  COMPUTATIONAL  SPACE 

m=l ; 

if  "sing  f  ield=  InitializeFields  (it ,  jt  ,kt)  ;  7.  SET-UP  A  CHUNK  OF  MEMORY 

else  field=  InitializeFieldsSingle(it , jt ,kt) ;  end 
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[C,D]=  InitializeMainGrid(C,D ,  it ,  jt  ,kt  ,m)  ;  7.  MAKE  ALL  CELLS  FREE  SPACE 

Is=  l:it; 

Js=  l:jt; 

Ks=  pb+l:kt; 
m=  2 ; 

[C,D]=  GetMaterial(C,D,  Is ,  Js  ,Ks,m)  ;  */.  MAKE  THE  CELLS  (Is,Js,Ks)  SUBSTRATE 

m=  1 ; 

[C,D] =  FillPML(C ,D ,dx ,dy , dz,pml ,pb , c , e , eO ,u,dt , it , jt ,kt ,m) ; 

7.  MAKE  ALL  OF  THE  PML  MATCHED  TO  FREE  SPACE 

m=  2 ; 

[C,D] =  changePML(C ,D,dx ,dy ,dz ,pml ,pb , c , e , eO ,u,dt , it , jt ,kt ,m,Ks (1) ) ; 

7.  MAKE  Ks  PML  MATCHED  TO  SUBSTRATE 

if  curve  7.  SPECIFY  SOURCE  POINT  AND  DISPLAY  VECTOR 

ks=  pb+round(kstruct/2) ; 
kd=  ks ; 

[field, is , j  s , id, jd, ANGLES] =  InitializePEC(f ield, it , jt ,kt ,pb,pml ,r ,dx,dy ,dz , Y,ks) ; 
7.  [field,  is ,  js ,  id,  jd,  ANGLES]  =  InitializePEC90  (field,  it ,  jt  ,kt  ,pb,pml  ,r  ,dx,dy  ,dz ,  Y,ks)  ; 

else 

7.  Location  of  current  source 
is=  pml+5; 
j  s=  pml+5 ; 
ks=  pb+3; 

7.  vector  to  display 
jd=  js+1; 
kd=  ks ; 

end 

fprintf  ( ’\n\nFDTD3D  preprocessing  took  7, .  2f  seconds  .  \n’  ,  toe)  ; 

7.  END  OF  PRE-PROCESSING— ONLY  DO  THE  ABOVE  ONCE 
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if  steps 


7.  BEGIN  TIME  STEPPING  LOOP 


last=  toe; 
flag=  1; 

for  t=l :TimeSteps 
if  flag 

if  t  >  TimeSteps/4  7  SEND  A  MESSAGE  WHEN  1/4  DONE 

quartertime=  toc-last; 

save ( [cd, ’ \Data\ ’ .name , ’quart ’] ,  ’ quart ertime ’ ) ; 
flag=  0; 

end; 

end; 

if  mod(t , round ( 1/ (f*dt)  )  )==0  7  SAVE  DATA  ONCE  A  PERIOD 

fprintf  (  ’  \nPeriod  7.. Of  of  7. Of  took  7.2f  minutes.  \n\t(7s)\n’ ,  .  .  . 

t /round (1/ (f *dt) ) , periods , (toc-last) /60, datestr (now) ) ; 
last=  toe; 
if  curve 

for  nn=  l:length(id) 

EZDATA(nn)=  field. Ez(id(nn) ,  jd(nn) ,kd) ; 

end; 

else  EZDATA(t/round(l/(f *dt) ) , : )=  field. Ez( : ,jd,kd) ;  end; 

end 

f ield=  UpdateE(f ield,C,D, it , jt ,kt ,DX,DY,DZ) ;  7  UPDATE  E  FIELDS 

field. Ez(is, js,ks)=  source(t);  7HARD  SOURCE 

if  curve  7  ZERO  THE  E  FIELDS  THIS  WAY  TO  REPRESENT  PEC 

field. Ex=  field. Ex  . *f ield.pecEx; 
field. Ey=  field. Ey  .  *f ield.pecEy ; 
field. Ez=  field. Ez  . *f ield.pecEz ; 

else  7  ZERO  THE  E  FIELDS  THIS  WAY 

field=  addPECTH (field,  it,  jt,  kt ,  pb,  pml) ; 

7  field=  addPECTF (field,  it,  jt,  kt ,  pb,  pml); 

7  f ield=  addPECMENZ (field,  it,  jt,  kt,  pb,  pml); 

end 
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f ield=  UpdateH(f ield,C,D, it , jt ,kt ,DX,DY,DZ) ; 


7  UPDATE  H  FIELDS 


if  mov  7  SAVE  A  FRAME  OF  FIELD  DATA 

f rame=  20;  7  SAVE  EVERY  frameTH  TIME  STEP 

figure (2) ; 
if  mod(t ,frame)==0 ; 

timestep=int2str (t) ; 

tview=squeeze (  field. Ez( : , : ,kd)  ); 

subplot ( ’position’ ,  [0 . 15  0.57  0.7  0.35]) 

imagesc(tview’ ) ; 

caxis  (  [-0 .1  0.1]); 

colorbar ; 

axis  image; 

title ([’Ez:  View  from  above  (+z) ;  time  step  =’ ,timestep] ) ; 

xlabel(’i  coordinate’); 
ylabelC’j  coordinate’); 

sview=squeeze (  field. Ez( : ,jd, : )  ); 

subplot ( ’position’ , [0 . 15  0.08  0.7  0.35]) 

imagesc (sview’ ) ; 

caxis  (  [-0 .1  0.1]); 

colorbar; 

axis  image; 

title([’Ez:  Cross-section  (-y) ;  time  step  =  ’ .timestep] ) ; 
xlabel(’i  coordinate’); 
ylabel(’k  coordinate’); 

tt=  t/frame; 

M(tt)=  getframe(gcf) ;  7  SAVE  A  FRAME  FOR  MOVIE 

if  av  7  SAVE  A  FRAME  FOR  . avi  MOVIE 

aviobj  =  addframe (aviobj ,M(tt) ) ; 

end 

end; 

end; 

end 
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if  mov 


movie(gcf ,M, 0 ,20) ;  "/„  Matlab  movie 

if  av  aviobj  =  close (aviobj ) ;  end  '/.  .avi  MOVIE  IS  NOW  IN  CURRENT  DIRECTORY 
save ( [cd , ’ \Data\M ’ , name] , ’ M ’ ) ; 
clear  M; 

end 

clear  field  /  DON’T  SAVE  ALL  OF  THE  FIELDS  OR  THE  .mat  FILE  IS  HUGE 

clear  C  D  source; 

End=  datestr(now) ; 

fprintf  (  ’  \nEnd  =  "/,s  .  \n\n’ ,End)  ; 

time=  toe; 

save( [cd, ’ \Data\ ’ ,name] ) ; 

end 
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InitializeFields.m 


function  field=  InitializeFields (it , jt ,kt) 

•/.  FUNCTION: 

•/.  1.  SET-UP  MEMORY  FOR  NUMBER  OF  CELLS:  (it.jt.kt)  OF  SIZE  double 

•/.  2.  CREATE  STRUCT  ’field’  THAT  CONTAINS  ALL  OF  THE  FIELD  COMPONENTS 


field. Ex=  zeros(it , jt+1 ,kt+l) ; 
field. Ey=  zeros(it+l , jt ,kt+l) ; 
field. Ez=  zeros (it+1 , jt+1 ,kt) ; 

field. Dx  =  zeros (it , jt+1 ,kt+l) ; 
field. Dy  =  zeros(it+l , jt,kt+l) ; 
field. Dz  =  zeros (it+1 , jt+1 ,kt) ; 
field. Dxn=  zeros (it , jt+1 ,kt+l) ; 
field. Dyn=  zeros(it+l , jt,kt+l) ; 
field. Dzn=  zeros (it+1 , jt+1 ,kt) ; 

field. Hx=  zeros(it+l , jt ,kt) ; 
field. Hy=  zeros(it , jt+1 ,kt) ; 
field. Hz=  zeros(it , jt ,kt+l) ; 


field. Bx  = 
field. By  = 
field. Bz  = 
field. Bxn= 
field. Byn= 
field. Bzn= 


zeros (it+1 , jt,kt) 
zeros (it , jt+1 ,kt) 
zeros (it , jt ,kt+l) 
zeros (it+1 , jt,kt) 
zeros (it , jt+1 ,kt) 
zeros (it , jt ,kt+l) 
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InitializeFieldsSingle.m 


function  field=  InitializeFields (it , jt ,kt) 

•/.  FUNCTION: 

•/.  1.  SET-UP  MEMORY  FOR  NUMBER  OF  CELLS:  (it.jt.kt)  OF  SIZE  single 

•/.  2.  CREATE  STRUCT  ’field’  THAT  CONTAINS  ALL  OF  THE  FIELD  COMPONENTS 


f ield .Ex=  single (zeros (it ,jt+l,kt+l)) ; 
f ield .Ey=  single (zeros (it +1 , jt ,kt+l) ) ; 
f ield .Ez=  single (zeros (it +1 , jt+1 ,kt) ) ; 

field. Dx  =  single (zeros (it , jt+1 ,kt+l) ) ; 
field. Dy  =  single (zeros (it+1 , jt ,kt+l) ) ; 
field. Dz  =  single (zeros (it+1 , jt+1 ,kt) ) ; 
f ield .Dxn=  single (zeros (it , jt+1 ,kt+l) ) ; 
f ield .Dyn=  single (zeros (it+1 , jt ,kt+l) ) ; 
f ield .Dzn=  single (zeros (it+1 , jt+1 ,kt) ) ; 

f ield.Hx=  single (zeros (it+1 , jt ,kt) ) ; 
f ield.Hy=  single (zeros (it , jt+1 ,kt) ) ; 
f ield .Hz=  single (zeros (it , jt ,kt+l) ) ; 


field. Bx  = 
field. By  = 
field. Bz  = 
field. Bxn= 
field.Byn= 
field.Bzn= 


single (zeros (it+1 , jt ,kt) ) 
single (zeros (it , jt+1 ,kt) ) 
single (zeros (it , jt ,kt+l) ) 
single (zeros (it+1 , jt ,kt) ) 
single (zeros (it , jt+1 ,kt) ) 
single (zeros (it , jt ,kt+l) ) 
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InitializeMainGrid.m 


function  [C,D]=  InitializeMainGrid(C ,D, it , jt ,kt ,m) 

•/.  FUNCTION: 

•/.  1.  SET-UP  ALL  CELLS  TO  CONTAIN  MATERIAL  ’m’ 


Ex=  zeros (it , jt+1 ,kt+l) 
Ey=  zeros (it+1 , jt ,kt+l) 
Ez=  zeros (it+1 , jt+1 ,kt) 
Hx=  zeros (it+1 , jt ,kt) ; 
Hy=  zeros (it , jt+1 ,kt) ; 
Hz=  zeros (it , jt ,kt+l) ; 


C . clEx= 

ones (size (Ex) ) 

C . c2Ex= 

ones (size (Ex) ) 

C . c3Ex= 

ones (size (Ex) ) 

C . c4Ex= 

ones (size (Ex) ) 

C . c5Ex= 

ones (size (Ex) ) 

C . c6Ex= 

ones (size (Ex) ) 

C . clEy= 

ones (size (Ey) ) 

C . c2Ey= 

ones (size (Ey) ) 

C . c3Ey= 

ones (size (Ey) ) 

C . c4Ey= 

ones (size (Ey) ) 

C . c5Ey= 

ones (size (Ey) ) 

C . c6Ey= 

ones (size (Ey) ) 

C . clEz= 

ones (size (Ez) ) 

C . c2Ez= 

ones (size (Ez) ) 

C . c3Ez= 

ones (size (Ez) ) 

C . c4Ez= 

ones (size (Ez) ) 

C . c5Ez= 

ones (size (Ez) ) 

C . c6Ez= 

ones (size (Ez) ) 

D . dlHx= 

ones (size (Hx) ) 

D . d2Hx= 

ones (size (Hx) ) 

D . d3Hx= 

ones (size (Hx) ) 

*C . cl (m) 
*C . c2(m) 
*C . c3(m) 
*C . c4(m) 
*C . c5(m) 
*C. c6(m) 

*C . cl (m) 
*C . c2(m) 
*C . c3(m) 
*C . c4(m) 
*C . c5(m) 
*C . c6(m) 

*C . cl (m) 
*C . c2(m) 
*C . c3(m) 
*C . c4(m) 
*C . c5(m) 
*C . c6(m) 

*D .  dl  (m) 
*D.d2(m) 
*D.d3(m) 
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D . d4Hx= 
D . d5Hx= 
D . d6Hx= 

D . dlHy= 
D . d2Hy= 
D . d3Hy= 
D . d4Hy= 
D . d5Hy= 
D . d6Hy= 

D . dlHz= 
D . d2Hz= 
D . d3Hz= 
D . d4Hz= 
D . d5Hz= 
D . d6Hz= 


ones (size (Hx) ) . *D . d4(m) 
ones (size (Hx) ) . *D . d5(m) 
ones (size (Hx) ) . *D . d6(m) 

ones (size (Hy) ) . *D . dl (m) 
ones (size (Hy) ) . *D . d2(m) 
ones (size (Hy) ) . *D . d3(m) 
ones (size (Hy) ) . *D . d4(m) 
ones (size (Hy) ) . *D . d5(m) 
ones (size (Hy) ) . *D . d6(m) 

ones (size (Hz) ) . *D . dl (m) 
ones (size (Hz) ) . *D . d2(m) 
ones (size (Hz) ) . *D . d3(m) 
ones (size (Hz) ) . *D . d4(m) 
ones (size (Hz) ) . *D . d5(m) 
ones (size (Hz) ) . *D . d6(m) 
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GetMaterial.m 


function  [C,D]=  GetMaterial (C ,D , I , J ,K ,mat) 


•/.  FUNCTION: 

•/.  1.  UPDATE  ’C’  AND  ’D’  TO  CHANGE 

C .  clEx(I ,  [J  max(J)+l],[K  max(K)+l])=  C. 
C.c2Ex(I,[J  max(J)+l],[K  max(K)+l])=  C. 
C .  c3Ex(I ,  [J  max(J)+l],[K  max(K)+l])=  C. 
C.c4Ex(I,[J  max(J)+l],[K  max(K)+l])=  C. 
C.c5Ex(I,[J  max(J)+l],[K  max(K)+l])=  C. 
C.c6Ex(I,[J  max(J)+l],[K  max(K)+l])=  C. 

C .  clEy (  [I  max(I)+l]  ,  J,  [K  max(K)+l])=  C. 
C .  c2Ey (  [I  max(I)+l]  ,  J,  [K  max(K)+l])=  C. 
C .  c3Ey (  [I  max(I)+l]  ,  J,  [K  max(K)+l])=  C. 
C .  c4Ey  (  [I  max(I)+l]  ,  J,  [K  max(K)+l])=  C. 
C .  c5Ey (  [I  max(I)+l]  ,  J,  [K  max(K)+l])=  C. 
C .  c6Ey (  [I  max(I)+l]  ,  J,  [K  max(K)+l])=  C. 

C .  clEz(  [I  max(I)+l],[J  max(J)+l]  ,K)=  C. 
C .  c2Ez(  [I  max(I)+l]  ,  [J  max(J)+l]  ,K)=  C. 
C .  c3Ez(  [I  max(I)+l],[J  max(J)+l]  ,K)=  C. 
C .  c4Ez(  [I  max(I)+l]  ,  [J  max(J)+l]  ,K)=  C. 
C .  c5Ez(  [I  max(I)+l]  ,  [J  max(J)+l]  ,K)=  C. 

C .  c6Ez(  [I  max(I)+l]  ,  [J  max(J)+l]  ,K)=  C. 

D .  dlHx(  [I  max(I)+l] , J,K)=  D.dl(mat); 
D.d2Hx([I  max(I)+l] , J,K)=  D.d2(mat); 
D.d3Hx([I  max(I)+l]  ,  J  ,K)=  D.d3(mat); 
D.d4Hx([I  max(I)+l] , J,K)=  D.d4(mat); 
D.d5Hx([I  max(I)+l] , J,K)=  D.d5(mat); 
D.d6Hx([I  max(I)+l]  ,  J  ,K)=  D.d6(mat); 


THE  MATERIAL  IN  CELLS  (I,J,K)  TO  ’mat’ 

cl (mat) ; 
c2(mat) ; 
c3(mat) ; 
c4(mat) ; 
c5(mat) ; 
c6(mat) ; 

cl (mat) ; 
c2(mat) ; 
c3(mat) ; 
c4(mat) ; 
c5(mat) ; 
c6(mat) ; 

cl (mat) ; 
c2(mat) ; 
c3(mat) ; 
c4(mat) ; 
c5(mat) ; 
c6(mat) ; 


D . dlHy (I , [J  max(J)+l] ,K)=  D.dl(mat); 
D.d2Hy(I,[J  max(J)+l] ,K)=  D.d2(mat); 
D.d3Hy(I,[J  max(J)+l] ,K)=  D.d3(mat); 
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D.d4Hy(I, [J 
D.d5Hy(I, [J 
D.d6Hy(I, [J 

D.dlHzU,  J, 
D.d2Hz(I, J, 
D.d3Hz(I, J, 
D.d4Hz(I, J, 
D.d5Hz(I, J, 
D.d6Hz(I, J, 


max(J)+l]  ,K)  = 
max(J)+l]  ,K)  = 
max(J)+l]  ,K)  = 

[K  max(K)  +  l]  )  = 
[K  max(K)  +  l]  )  = 
[K  max(K)  +  l]  )  = 
[K  max(K)  +  l]  )  = 
[K  max(K)  +  l]  )  = 
[K  max(K)  +  l]  )  = 


D.d4(mat) ; 
D.d5 (mat) ; 
D.d6 (mat) ; 

D . dl (mat) ; 
D.d2 (mat) ; 
D.d3(mat) ; 
D . d4(mat) ; 
D.d5 (mat) ; 
D.d6 (mat) ; 
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FillPML.m 


function  [C ,D] =  FillPML(C,D,dx,dy ,dz ,pml ,pb, c , e ,eO,u,dt , it , jt ,kt ,m) 

'/,  Program  author:  Keely  J.  Willis,  Ph.D.  Student 

'/,  UW  Computational  Electromagnetics  Laboratory 

'/,  Director:  Prof.  Susan  C.  Hagness 

•/.  FUNCTION: 

•/.  1.  FILL  THE  OUTSIDE  ’pml’  NUMBER  OF  CELLS  AS  A  POLYNOMIAL  GRADED  UNIAXIAL 

•/.  PERFECTLY  MATCHED  LAYER  (UPML)  MATCHED  TO  FREE  SPACE  BY  UPDATING  ’C’  AND  ’D’ 

7 

rmax=  exp(-16);  7.  MAX  DESIRED  NORMAL  REFLECTION 

order=  4;  7  ORDER  OF  GRADING 

'/,  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-varying  material  properties 

d=  pml*dx;  7.  <x>  THICKNESS  OF  PML  (m) 

sigmaMax=  - (order+1) *log(rmax) *eO*c* . 5/d;  7  MAX  <x>  CONDUCTIVITY 

sigf actor=sigmaMax/ (dx* (d~order) * (order+1) ) ; 

for  i=l:pml 

xl=(pml-i+l)*dx; 

x2=(pml-i)*dx; 

sigma=sigf actor* (xl~ (order* l)-x2~ (order+1) ) ; 
facm=(2*e(m)-sigma*dt) ; 
facp=(2*e(m)+sigma*dt) ; 

7  MATERIAL  COEFFICIENTS  FOR  CELLS  IN  THE  CENTER  OF  THE  UPML 
C . c5Ex(i , : , : )=f acp; 

C . c5Ex(it-i+l , : , : )=f acp ; 

C . c6Ex(i , : , : )=f acm; 

C .  c6Ex(it-i+l , : , : )=f acm; 

D .  dlHz (i , : , : )=f acm/f acp ; 

D . dlHz (it-i+1 , : , : )=f acm/f acp ; 

D . d2Hz(i , : , : )=2*e (m) *dt/f acp ; 

D . d2Hz (it-i+1 , : , : )=2*e (m) *dt/f acp; 
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D . d3Hy ( i , : , : ) =f  acm/ f  acp ; 

D . d3Hy (it-i+1 , : , : )=f acm/f acp ; 

D . d4Hy ( i , : , : )  =  1 /f  acp/ u (m) ; 

D . d4Hy (it-i+1 , : , : )=l/f acp/u(m) ; 

xl=(pml-i+l . 5) *dx ; 
x2=(pml-i+0 . 5) *dx ; 

sigma=sigf actor* (xl~ (order* l)-x2~ (order* 1) ) ; 
facm=(2*e(m)-sigma*dt) ; 
facp=(2*e(m)+sigma*dt) ; 

"/.  MATERIAL  COEFFICIENTS  FOR  CELLS  ON  THE  UPML  BOUNDARY 
C . clEz (i , : , : )=f acm/f acp ; 

C . clEz (it+l-i+1 , : , : )=f acm/f acp ; 

C . c2Ez (i , : , : )=2*e (m) *dt/f acp ; 

C . c2Ez (it+l-i+1 , : , : )=2*e(m) *dt/f acp; 

C . c3Ey ( i , : , : ) =f acm/ f  acp ; 

C . c3Ey (it+l-i+1 , : , : )=f acm/f acp ; 

C . c4Ey ( i , : , : )  =  1 /f  acp/ e (m) ; 

C .  c4Ey (it+l-i+1 , : , : )=l/f acp/e(m) ; 

D .  d5Hx ( i , : , : ) =f acp ; 

D . d5Hx( it+l-i+1 , : , : )=f acp ; 

D . d6Hx ( i , : , : ) =f acm ; 

D . d6Hx( it+l-i+1 , : , : )=f acm; 

end 


'/•  yyyyyyyyyyyyyyyyyyyyyyyyyyyyy_varying  material  properties 

d=  pml*dy ;  7,  <y>  THICKNESS  OF  PML  (m) 

sigmaMax=  -  (order+1)  *log(rmax)  *eO*c*  .  5/d;  "/,  MAX  <y>  CONDUCTIVITY 

sigf actor=  sigmaMax/ (dy* (d~order) * (order+1) ) ; 

for  j=l:pml 

yl=  (pml-j+l)*dy; 
y2=  (pml-j ) *dy ; 

sigma=  sigf actor* (yl~ (order+1) -y2~ (order+1) ) ; 
facm=  (2*e (m) -sigma*dt) ; 
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facp=  (2*e(m)+sigma*dt) ; 


"/.  MATERIAL  COEFFICIENTS  FOR  CELLS  IN  THE  CENTER  OF  THE  UPML 
C .  c5Ey  (  :  ,  j  ,  :  )  =f  acp ; 

C . c5Ey ( : , jt-j+1 , : )=f acp ; 

C .  c6Ey  (  :  ,  j  ,  :  )  =f  acm ; 

C .  c6Ey ( : , jt-j+1 , : )=f acm; 

D .  dlHx ( : , j , : ) =f acm/ f  acp ; 

D . dlHx( : , jt-j+1 , : )=f acm/f acp ; 

D . d2Hx( : , j , : )=2*e (m) *dt/f acp ; 

D . d2Hx( : , jt-j+1 , : )=2*e(m) *dt/f acp; 

D . d3Hz ( : , j , : ) =f acm/ f  acp ; 

D . d3Hz( : , jt-j+1 , : )=f acm/f acp ; 

D . d4Hz ( : , j , : ) =l/f acp/ u (m) ; 

D . d4Hz( : , jt-j+1 , : )=l/f acp/u(m) ; 

yl=  (pml-j+1 .5)*dy; 
y2=  (pml-j+0.5)*dy; 

sigma=  sigf actor* (yl~ (order+1) -y2~ (order+1) ) ; 
facm=  (2*e (m) -sigma*dt) ; 
facp=  (2*e(m)+sigma*dt) ; 

"/.  MATERIAL  COEFFICIENTS  FOR  CELLS  ON  THE  UPML  BOUNDARY 
C . clEx( : , j , : )=f acm/f acp ; 

C . clEx( : , jt+l-j+1 , : )=f acm/f acp; 

C . c2Ex( : , j , : )=2*e (m) *dt/f acp ; 

C . c2Ex( : , jt+l-j+1, : )=2*e(m)*dt/f acp ; 

C . c3Ez ( : , j , : ) =f acm/ f  acp ; 

C . c3Ez (: , jt+l-j+1, : )=f acm/f acp ; 

C . c4Ez( : , j , : )=l/f acp/e(m) ; 

C .  c4Ez (: , jt+l-j+1, : )=l/f acp/e(m) ; 

D .  d5Hy ( : , j , : ) =f acp ; 

D . d5Hy ( : , jt+l-j+1, : )=f acp ; 

D . d6Hy ( : , j , : ) =f  acm ; 

D . d6Hy ( : , jt+l-j+1, : )=f acm; 

end 
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'/,  zzzzzzzzzzzzzzzzzzzzz-varying  material  properties 

d=  pml*dz;  */.  <z>  THICKNESS  OF  PML  (m) 

sigmaMax=  - (order+1) *log(rmax) *eO*c* . 5/d;  "/,  MAX  <z>  CONDUCTIVITY 

sigf actor=  sigmaMax/ (dy* (d~order) * (order+1) ) ; 

for  k=l:pml 

zl=(pml-k+l)*dz; 

z2=(pml-k)*dz; 

sigma=sigf actor* (zl~ (order* l)-z2~ (order+1) ) ; 
facm=(2*e(m)-sigma*dt) ; 
facp=(2*e(m)+sigma*dt) ; 

*/.  MATERIAL  COEFFICIENTS  FOR  CELLS  IN  THE  CENTER  OF  THE  UPML 


C . c5Ez( : 

, : ,k)=facp; 

C . c6Ez( : 

, : ,k)=facm; 

D . dlHy ( : 

, : ,k)=f acm/f acp; 

D.d2Hy(: 

, : ,k)=2*e(m)*dt/facp; 

D.d3Hx(: 

, : ,k)=f acm/f acp; 

D.d4Hx(: 

, : ,k)=l/f acp/u(m) ; 

zl=(pml-k+l . 5) *dz ; 
z2=(pml-k+0 . 5) *dz ; 

sigma=sigf actor* (zl~ (order* l)-z2~ (order+1) ) ; 
facm=(2*e(m)-sigma*dt) ; 
facp=(2*e(m)+sigma*dt) ; 

"/.  MATERIAL  COEFFICIENTS  FOR  CELLS  ON  THE  UPML  BOUNDARY 
C . clEy ( : , : ,k)=f acm/f acp ; 

C . c2Ey ( : , : ,k)=2*e (m) *dt/f acp ; 

C . c3Ex ( : , : , k) =f  acm/ f  acp ; 

C . c4Ex( : , : ,k)=l/f acp/e(m) ; 

D  .  d5Hz  ( :  ,  :  ,  k)  =f  acp ; 

D  .  d6Hz  ( :  ,  :  ,  k)  =f  acm ; 

end 
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y.  <x>  NORMAL  PEC  WALLS 
C.clEz(l, : , : )=-l . 0 ; 

C . clEz(it+l , : , : )=-l . 0 ; 

C . c2Ez(l , : , : )=0.0; 
C.c2Ez(it+l, : , :)=0.0; 

C . c3Ey (1 , : , :)=-1.0; 

C . c3Ey (it+1 , : , : )=-l . 0 ; 

C . c4Ey (1 , : , : )=0.0; 

C . c4Ey (it+1 , : , :)=0.0; 

•/.  <y>  NORMAL  PEC  WALLS 


C.clEx(: 

C.clEx(: 

,jt+l, :)=-l; 

C . c2Ex( : 

,1,0=0; 

C . c2Ex( : 

,jt+l,  0=0; 

C . c3Ez( : 

,1,0=-1; 

C . c3Ez( : 

,jt+l,  0=-l; 

C . c4Ez  (  : 

,1, 0=0; 

C . c4Ez ( : 

,  jt+1,  0=0; 

•/.  <z> 

NORMAL  PEC  WALLS 

C.clEy(: 

> : > l)=-l ; 

C . c2Ey ( : 

, : ,1)=0; 

C . c3Ex  (  : 

>  :  >  l)  =— l ; 

C . c4Ex ( : 

, : ,1)=0; 
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ChangePML.m 


function  [C ,D] =  changePML(C,D ,dx,dy ,  dz ,pml ,pb, c ,e , eO ,u,dt , it , jt ,kt ,m,K) 

'/,  Program  author:  Keely  J.  Willis,  Ph.D.  Student 

'/,  UW  Computational  Electromagnetics  Laboratory 

'/,  Director:  Prof.  Susan  C.  Hagness 

•/.  FUNCTION: 

•/.  1.  CHANGE  PML  OF  OUTSIDE  CELLS  OF  LAYERS  (K:kt)  TO  BE  MATCHED  TO  MATERIAL 

•/.  ’m’  BY  UPDATING  ’C’  AND  >D’ 

7 

rmax=  exp(-16);  7.  MAX  DESIRED  NORMAL  REFLECTION 

order=  4;  7,  ORDER  OF  GRADING 

'/,  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-varying  material  properties 
d=  pml*dx;  7.  <x>  THICKNESS  OF  PML  (m) 

sigmaMax=  - (order+1) *log(rmax) *eO*c* . 5/d;  "/,  MAX  <x>  CONDUCTIVITY 

sigf actor=sigmaMax/ (dx* (d~order) * (order+1) ) ; 

for  i=l:pml 

xl=  (pml-i+l)*dx; 
x2=  (pml-i)*dx; 

sigma=  sigf actor* (xl“ (order+1) -x2“ (order+1) ) ; 
f acm=  (2*e0-sigma*dt) ; 
f acp=  (2*e0+sigma*dt) ; 

7  MATERIAL  COEFFICIENTS  FOR  CELLS  IN  THE  CENTER  OF  THE  UPML 
C . c5Ex ( i , : , K : end) =f acp ; 

C . c5Ex(it-i+l , : ,K:end)=facp; 

C . c6Ex ( i , : , K : end) =f acm ; 

C .  c6Ex(it-i+l , : ,K:end)=facm; 

D. dlHz(i, : ,K:end)=facm/facp; 

D . dlHz (it-i+1 , : ,K:end)=facm/facp; 

D.d2Hz(i, : ,K:end)=2*e(m)*dt/facp; 

D . d2Hz(it-i+l , : ,K:end)=2*e(m)*dt/facp; 
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D . d3Hy ( i , : , K : end) =f acm/f  acp ; 

D . d3Hy (it-i+1 , : ,K:end)=f acm/f acp; 

D . d4Hy (i , : ,K : end)=l/f acp/u(m) ; 

D . d4Hy (it-i+1 , : ,K : end)=l/f acp/u(m) ; 

xl=  (pml-i+1 .5)*dx; 
x2=  (pml-i+0.5)*dx; 

sigma=  sigf actor* (xl~ (order+1) -x2~ (order+1) ) ; 
f acm=  (2*e0-sigma*dt) ; 
f acp=  (2*e0+sigma*dt) ; 

"/.  MATERIAL  COEFFICIENTS  FOR  CELLS  ON  THE  UPML  BOUNDARY 
C.clEz(i, : ,K: end) =  f  acm/ f  acp ; 

C . clEz (it+l-i+1 , : ,K:end)=  facm/facp; 

C.c2Ez(i, : ,K:end)=  2*e(m)*dt/facp; 

C . c2Ez (it+l-i+1 , : ,K:end)=  2*e(m) *dt/facp; 

C . c3Ey ( i , : , K : end) =  f  acm/ f  acp ; 

C . c3Ey (it+l-i+1 , : ,K:end)=  facm/facp; 

C . c4Ey ( i , : , K : end) =  1 /f  acp/e (m) ; 

C .  c4Ey (it+l-i+1 , : ,K : end)=  1/f acp/ e (m) ; 

D. d5Hx(i,:,K:end)=  facp; 

D . d5Hx(it+l-i+l , : ,K : end)=  facp; 

D . d6Hx ( i , : , K : end) =  f  acm ; 

D . d6Hx(it+l-i+l , : ,K : end)=  facm; 

end 


'/•  yyyyyyyyyyyyyyyyyyyyyyyyyyyyy_varying  material  properties 

d=  pml*dy ;  */.  <y>  THICKNESS  OF  PML  (m) 

sigmaMax=  -  (order+1)  *log(rmax)  *eO*c*  .  5/d;  '/,  MAX  <y>  CONDUCTIVITY 

sigf actor=  sigmaMax/ (dy* (d~order) * (order+1) ) ; 

for  j=l:pml 

yl=  (pml-j+l)*dy; 
y2=  (pml-j ) *dy ; 

sigma=  sigf actor* (yl~ (order+1) -y2~ (order+1) ) ; 
f acm=  (2*e0-sigma*dt) ; 
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f acp=  (2*e0+sigma*dt) ; 


"/.  MATERIAL  COEFFICIENTS  FOR  CELLS  IN  THE  CENTER  OF  THE  UPML 
C . c5Ey ( : , j ,K : end)=  facp; 

C . c5Ey ( : , jt- j+1 ,K : end)=  facp; 

C . c6Ey ( : , j ,K : end)=  facm; 

C .  c6Ey ( : , jt- j+1 ,K : end)=  facm; 

D .  dlHx( : , j ,K : end)=  facm/facp; 

D . dlHx( : , jt- j+1 ,K : end)=  facm/facp; 

D . d2Hx( : , j ,K : end)=  2*e(m) *dt/f acp; 

D . d2Hx( : , jt- j+1 ,K : end)=  2*e (m) *dt/f acp; 

D . d3Hz( : , j ,K : end)=  facm/facp; 

D . d3Hz( : , jt- j+1 ,K : end)  =  facm/facp; 

D . d4Hz( : , j ,K : end)=  l/facp/u(m) ; 

D . d4Hz( : , jt- j+1 ,K : end)=  1/f acp/u(m) ; 

yl=  (pml-j+1 .5)*dy; 
y2=  (pml-j+0.5)*dy; 

sigma=  sigf actor* (yl~ (order+1) -y2~ (order+1) ) ; 
f acm=  (2*e0-sigma*dt) ; 
f acp=  (2*e0+sigma*dt) ; 

"/.  MATERIAL  COEFFICIENTS  FOR  CELLS  ON  THE  UPML  BOUNDARY 
C . clEx( : , j ,K : end)=  facm/facp; 

C . clEx( : , jt+l-j+1 ,K: end)=  facm/facp ; 

C . c2Ex( : , j ,K : end)=  2*e(m) *dt/f acp; 

C . c2Ex( : , jt+l-j+1, K: end)=  2*e(m) *dt/f acp; 

C . c3Ez ( : , j ,K : end)=  facm/facp; 

C . c3Ez (: , jt+l-j+1, K: end)=  facm/facp ; 

C . c4Ez( : , j ,K : end)=  l/facp/e(m); 

C .  c4Ez (: , jt+l-j+1, K: end)=  1/f acp/ e (m) ; 

D .  d5Hy ( : , j , K : end) =  f  acp ; 

D . d5Hy( :, jt+l-j+1 ,K : end)=  facp; 

D . d6Hy ( : , j ,K : end)=  facm; 

D . d6Hy( :, jt+l-j+1 ,K : end)=  facm; 

end 
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InitializePEC.m 


function  [field, is , js , id, jd, ANGLES] = 

InitializePEC(f ield, it , jt ,kt ,pb,pml ,r ,dx ,dy ,dz , Y,ks) 

7.  FUNCTION: 

•/.  1.  ASSIGN  LOCATION  OF  PEC  FOR  180  DEGREE  CURVE  BY  ZEROING  TANGENTIAL  E  FIELDS 

•/.  2.  SPECIFY  THE  SOURCE  CELL 

•/.  3.  SPECIFY  CELLS  TO  BE  RECORDED  FOR  PROPAGATION  CONSTATNT  EXTRACTION 

•/. 

•/.  OUTPUTS : 

•/.  ‘field. pecEx’  ,  ‘field. pecEy’  ,  ‘field. pecEz’  :  SPECIFY  CELLS  OF  E  FIELDS  TO  BE 

•/.  ZEROED  TO  MODEL  PEC 

•/. 

•/.  ‘ANGLES’:  A  VECTOR  OF  ANGLES  TO  EACH  (id,jd)  CELL  FROM  CENTER  OF  CIRCLE 

•/. 

•/.  ‘id’  AND  ‘jd’:  SPECIFY  THE  FIELD  DATA  TO  BE  SAVED 

•/. 

•/.  ‘is’  AND  ‘js’:  SPECIFY  THE  SOURCE  CELL(S) 

7. 

PEC=  0; 

•/.  CREATE  MEMORY  LOCATIONS  FOR  PEC 
f ield .pecEx=  single (ones (it ,jt+l,kt+l)) ; 
f ield .pecEy=  single (ones (it +1 , jt ,kt+l) ) ; 
f ield .pecEz=  single (ones (it +1 , jt+1 ,kt) ) ; 

'/,  Put  ground  plate  on  entire  floor  of  computational  space 
Igx=  l:it; 

Igy=  1: (it+1) ; 

Jgx=  1 : (jt+1) ; 

Jgy=  1 : jt ; 

Kg=  (kt+1) ; 

field. pecEx(Igx, Jgx,Kg)=  PEC; 
field. pecEy(Igy, Jgy,Kg)=  PEC; 
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Kc=  pb+1 ; 
Kw=  Kc:kt; 
INSIDE=  1; 


7.  vertical  location  of  strip 
7.  vertical  location  of  wall 
7.  select  inside  or  outside  wall 


'/,  (distOX,distOY)  are  the  coordinates  (m)  of  the  center  of  the  "circle" 
distOY=  r/sqrt(2)  +jt*dy  -dy*(2+pml); 
distOX=  round(it/2)*dx; 

id=  0;  jd=  0;  ANGLES=  0; 
for  x=  l:it 

for  y=  1 : jt 

'/,  (distX.distY)  are  the  coordinates  (m)  of  the  current  cell  at  ’  x’  and  ’y’ 
distX=  (x-l)*dx; 
distY=  (y-l)*dy; 

'/,  DIST  is  the  distance  from  center  of  the  "circle"  to  (distX,distY) 

DIST=  sqrt ( (distOX-distX) ~2  +(distOY-distY) ~2) ; 

if  (DIST  >  r)  &  (DIST  <  (r+Y) ) 
field. pecEx(x,y,Kc)=  PEC; 
field. pecEy(x,y,Kc)=  PEC; 

end; 

if  INSIDE  7.  PUT  WALL  ON  INSIDE  CURVE 

if  abs (DIST-r)  <  (dy/2) 

f ield.pecEx(x,y ,Kw)=  PEC; 
field. pecEz(x,y,Kw)=  PEC; 

end; 

else  7.  PUT  WALL  ON  OUTSIDE  CURVE 

if  abs (r+Y-DIST)  <  (dy/2) 

f ield.pecEx(x,y ,Kw)=  PEC; 
field. pecEz(x,y,Kw)=  PEC; 

end; 

end 

if  INSIDE  7.  PUT  SOURCE  CELL  AND  DISPLAY  VECTOR  NEAR  OUTSIDE  EDGE 


'/.Concucting  strip 
'/.Concucting  strip 
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is=  pml+4; 
js=  jt-pml-4; 

if  abs (DIST- (r+Y-2*dy) )  <  (dy/2) 

ANGLES=  [ANGLES  atan2(dist0Y-distY,dist0X-distX)*180/pi] ; 
id=  [id  x] ; 
jd=  [jd  y]  ; 

end; 


else  "/.  PUT  SOURCE  CELL  AND  DISPLAY  VECTOR  NEAR  INSIDE  EDGE 

is=  pml+Y-2; 
js=  jt-pml-4; 

if  abs (DIST- (r+2*dy) )  <  (dy/2) 

ANGLES=  [ANGLES  atan2(dist0Y-distY,dist0X-distX)*180/pi] ; 
id=  [id  x] ; 
jd=  [jd  y]  ; 

end; 

end 

end; 

end; 

id ( 1)  =  []  ;  jd(l)  =  []  ;  ANGLES ( 1)  =  []  ; 

if  id==  []  |  jd==  []  error  (’YOU  W0N”T  GET  DATA,  DUMMY’);  end; 
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InitializePEC90.m 


function  [field, is ,js , id, jd, ANGLES] =  ... 

InitializePEC90(f ield, it , jt ,kt ,pb,pml ,r ,dx,dy ,dz, Y,ks) 

7.  FUNCTION: 

•/.  1.  ASSIGN  LOCATION  OF  PEC  FOR  90  DEGREE  CURVE  BY  ZEROING  TANGENTIAL  E  FIELDS 

•/.  2.  SPECIFY  THE  SOURCE  CELL 

•/.  3.  SPECIFY  CELLS  TO  BE  RECORDED  FOR  PROPAGATION  CONSTATNT  EXTRACTION 

•/. 

•/.  OUTPUTS : 

•/.  ‘field. pecEx’  ,  ‘field. pecEy’  ,  ‘field. pecEz’  :  SPECIFY  CELLS  OF  E  FIELDS  TO  BE 

•/.  ZEROED  TO  MODEL  PEC 

•/. 

•/.  ‘ANGLES’:  A  VECTOR  OF  ANGLES  TO  EACH  (id,jd)  CELL  FROM  CENTER  OF  CIRCLE 

•/. 

•/.  ‘id’  AND  ‘jd’:  SPECIFY  THE  FIELD  DATA  TO  BE  SAVED 

•/. 

•/.  ‘is’  AND  ‘js’:  SPECIFY  THE  SOURCE  CELL(S) 

7. 

PEC=  0; 

•/.  CREATE  MEMORY  LOCATIONS  FOR  PEC 
field. pecEx=  ones (it , jt+1 ,kt+l) ; 
field. pecEy=  ones (it+1 , jt ,kt+l) ; 
field. pecEz=  ones (it+1 , jt+1 ,kt) ; 

'/,  Put  ground  plate  on  entire  floor  of  computational  space 
Igx=  1 : it ; 

Igy=  1: (it+1) ; 

Jgx=  1:  (jt+1) ; 

Jgy=  l:jt; 

Kg=  (kt+1) ; 

field. pecEx(Igx, Jgx,Kg)=  PEC; 
field. pecEy(Igy, Jgy,Kg)=  PEC; 
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Kc=  pb+1; 
Kw=  Kc:kt; 
INSIDE=  1; 


7.  vertical  location  of  strip 
7.  vertical  location  of  wall 
'/„  select  inside  or  outside  wall 


'/,  (distOX,distOY)  are  the  coordinates  (m)  of  the  center  of  the  "circle" 
distOY=  r/sqrt(2)  +jt*dy  -dy*(2+pml); 
distOX=  round(it/2)*dx; 

id=  0;  jd=  0;  ANGLES=  0; 
for  x=  1 : it 

for  y=  l:jt+l 

'/,  (distX.distY)  are  the  coordinates  (m)  of  the  current  cell  at  ’x’  and  ’y’ 
distX=  (x-l)*dx; 
distY=  (y-l)*dy; 

'/,  DIST  is  the  distance  from  center  of  the  "circle"  to  (distX,distY) 

DIST=  sqrt ( (distOX-distX) ~2  +(distOY-distY) ~2) ; 

if  (DIST  >  r)  &  (DIST  <  (r+Y) ) 
field. pecEx(x,y,Kc)=  PEC; 
field. pecEy(x,y,Kc)=  PEC; 

end; 

if  INSIDE  7.  PUT  WALL  ON  INSIDE  CURVE 

if  abs (DIST-r)  <  (dy/2) 

f ield.pecEx(x,y ,Kw)=  PEC; 
field. pecEz(x,y,Kw)=  PEC; 

end; 

else  7,  PUT  WALL  ON  OUTSIDE  CURVE 

if  abs (r+Y-DIST)  <  (dy/2) 

f ield.pecEx(x,y ,Kw)=  PEC; 
field. pecEz(x,y,Kw)=  PEC; 

end; 

end 

if  INSIDE  7  PUT  SOURCE  CELL  AND  DISPLAY  VECTOR  NEAR  OUTSIDE  EDGE 

if  x==(pml+5) 


'/.Concucting  strip 
'/.Concucting  strip 
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if  abs (DIST- (r+Y-2*dy) )  <  (dy/2) 
is=  x; 

js=  y; 

end; 

end; 

if  abs (DIST- (r+Y-2*dy) )  <  (dy/2) 

ANGLES=  [ANGLES  atan2(dist0Y-distY,dist0X-distX)*180/pi] ; 
id=  [id  x] ; 
jd=  [jd  y]  ; 

end; 

else  "/.  PUT  SOURCE  CELL  AND  DISPLAY  VECTOR  NEAR  INSIDE  EDGE 

if  x==(pml+5) 

if  abs (DIST- (r+Y-2*dy))  <  (dy/2) 
is=  x; 

js=  y; 

end; 

end; 

if  abs (DIST- (r+2*dy) )  <  (dy/2) 

ANGLES=  [ANGLES  atan2(dist0Y-distY,dist0X-distX)*180/pi] ; 
id=  [id  x] ; 
jd=  [jd  y]  ; 

end; 

end 

end; 

end; 

id(l)  =  []  ;  jd(l)  =  []  ;  ANGLES(1)  =  []; 

if  id==  []  |  jd==  []  error  (’YOU  W0N”T  GET  DATA,  DUMMY’);  end; 


119 


UpdateE.m 


function  field=  UpdateE(f ield.C ,D, it , jt ,kt ,DX ,DY,DZ) 

7. 

•/.  FUNCTION: 

•/.  1.  UPDATE  E  FIELD  COMPONENT  OF  ALL  CELLS  FOR  ONE  TIME  STEP  BASED  ON  EXAMPLE  OF 

7.  EQ.  7.74  OF  COMPUTAIONAL  ELECTRODYNAMICS  BY  TAFLOVE/HAGNESS  2ND  ED. 

•/. 


f ield .Dxn( :,2:jt,2:kt)=  C. clEx( : ,2:jt,2:kt) . *f ield. Dx( : ,2:jt,2:kt) . . . 

+C.c2Ex(: ,2:jt,2:kt) .*(  (field.Hz(: ,2 : jt ,2 :kt) -f ield. Hz( : , 1 : jt-1 ,2 :kt) )*DY  ... 

-(f ield.Hy ( : ,2 : jt , 2 :kt) -f ield. Hy ( : ,2 : jt , 1 :kt-l) )*DZ  ) ; 

f ield .Dyn(2 : it , : ,2 :kt)=  C . clEy (2 : it , : , 2 :kt) . *f ield.Dy (2 : it , : ,2 :kt) . . . 

+  C . c2Ey (2:it, : ,2:kt) .*(  (f ield.Hx(2 : it , : ,2 :kt)-f ield.Hx(2 : it , : , 1 :kt-l) ) *DZ  ... 

-(field. Hz (2: it, : ,2 :kt)-f ield. Hz (1 : it-1 , : ,2:kt))*DX  ) ; 

f ield.Dzn(2 : it ,2 : jt , : )=  C. clEz(2 : it , 2 : jt , : ) . *f ield.Dz (2 : it , 2 : jt , : ) . . . 

+  C . c2Ez (2 : it , 2 : jt , :) .*(  (f ield.Hy (2 : it , 2 : jt )-f ield.Hy (1 : it-1 ,2 : jt ,:)) *DX  ... 

- (f ield. Hx(2 : it ,2 : jt , :) -field. Hx(2 : it , 1 : jt-1 , :))*DY  ) ; 

field.Ex(:,2:jt,2:kt)=  C. c3Ex( : ,2:jt,2:kt) . *f ield. Ex ( : ,2:jt,2:kt) . . . 

+  C . c4Ex( :,2:jt,2:kt).*(  C. c5Ex( : ,2:jt,2:kt) . *f ield.Dxn( :,2:jt,2:kt)... 

-C . c6Ex( : ,2:jt,2:kt) . *f ield. Dx( : ,2 : jt ,2 :kt)  ); 

f ield .Ey (2 : it , : , 2 : kt)=  C . c3Ey (2 : it , : ,2 :kt) . *f ield.Ey (2 : it , : ,2 :kt) . . . 

+  C . c4Ey (2 : it , : , 2 :kt) . * (  C . c5Ey (2 : it , : ,2 :kt) . *f ield.Dyn(2 : it , : ,2 :kt) . . . 

-C . c6Ey (2 : it , : ,2 :kt) . *f ield.Dy (2 : it , : ,2 :kt)  ) ; 

field.Ez(2:it,2:jt,:)=  C. c3Ez(2 : it , 2 : jt , : ) . *f ield.Ez (2:it,2:jt, :) . . . 

+  C . c4Ez (2 :it,2:jt,:).*(  C. c5Ez (2 : it , 2 : jt , : ) . *f ield.Dzn(2 : it , 2 : jt , : ) . . . 

-C . c6Ez(2 : it , 2 : jt , : ) . *f ield.Dz (2:it,2:jt,:)  ); 

f ield. Dx=f ield. Dxn; 
f ield. Dy=f ield. Dyn; 
f ield. Dz=f ield. Dzn; 
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UpdateH.m 


function  field=  UpdateH(f ield,C ,D, it , jt ,kt ,DX ,DY,DZ) 

7. 

•/.  FUNCTION: 

•/.  1.  UPDATE  H  FIELD  COMPONENT  OF  ALL  CELLS  FOR  ONE  TIME  STEP  BASED  ON  EXAMPLE  OF 

7.  EQ.  7.74  OF  COMPUTAIONAL  ELECTRODYNAMICS  BY  TAFLOVE/HAGNESS  2ND  ED. 

•/. 

f ield .Bxn(2 : it , : , : D . dlHx(2 : it *f ield.Bx (2 : it . 

-  D . d2Hx (2 : it (  (f ield.Ez(2 : it ,2 : jt+1 , : )-f ield. Ez(2 : it , 1 : jt , : ) ) *DY  ... 

- (f ield.Ey (2 : it , : ,2 :kt+l)-f ield.Ey (2: it , : ,l:kt))*DZ  ) ; 

f ield .Byn( : , 2 : jt , : )=  D.dlHyC : ,2: jt , : ) . *  f ield. By ( : , 2 : jt , : ) . . . 

-  D.d2Hy(: ,2: jt, :) .*(  (f ield. Ex ( : , 2 : jt , 2 :kt+l) -f ield. Ex( : , 2 : jt , 1 :kt) ) *DZ  ... 

- (f ield.Ez(2 : it+1 , 2 : jt , : )-f ield. Ez(l : it ,2: jt , :))*DX  ) ; 

f ield .Bzn( : , : ,2 :kt)=  D . dlHz ( : , : ,2 : kt) . *  f ield. Bz( : , : ,2 :kt) . . . 

-  D.d2Hz(: , : ,2:kt) .*(  (f ield.Ey (2 : it+1 , : ,2 : kt) -f ield. Ey (1 : it , : , 2 :kt) ) *DX  ... 

- (field. Ex( : , 2 : jt+1 , 2 :kt) -f ield. Ex( : , 1 : jt , 2 :kt) ) *DY  ) ; 

field. Hx(2:it, :  ,  :)  =  D.d3Hx(2:it, : , . *f ield. Hx(2 : it , : , D.d4Hx(2:it, 

. * (D . d5Hx(2 : it , : , : ) . *f ield. Bxn(2 : it , : , : )  -  D . d6Hx(2 : it , : , : ) . *f ield.Bx (2 : it , : , : ) ) ; 

field.HyC: ,2: jt, :)=  D.d3Hy(: ,2: jt, :) ,*field.Hy(: ,2: jt, :)+  D.d4Hy(: ,2: jt, :) . . . 

. * (D . d5Hy ( : , 2 : jt , : ) . *f ield. Byn( : ,2 : jt , : )  -  D. d6Hy ( : ,2:jt, :) ,*field.By(: ,2:jt, :)) ; 

field. Hz ( : , : ,2:kt)=  D.d3Hz(: , : ,2:kt) ,*field.Hz(: , : ,2:kt)+  D.d4Hz(: , : ,2:kt) . . . 

. *(D.d5Hz( : , : ,2:kt) . *f ield. Bzn( : , : ,2:kt)  -  D.d6Hz(: , : ,2:kt) ,*field.Bz(: , : ,2:kt)) ; 

f ield. Bx=f ield. Bxn; 
field . By=f ield . Byn ; 
f ield. Bz=f ield. Bzn; 
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AddPECTH.m 


function  field=  addPECTH (field,  it,  jt,  kt ,  pb,  pml) 

7. 

'/.  FUNCTION: 

'/.  1.  ASSIGN  ZERO  THICKNESS  PEC  FOR  A  THW  ANTENNA  TO  BY  ZEROING 

'/.  TANGENTIAL  E  FIELDS  CORRESPONDING  TO  CONDUCTOR  LOCATIONS 

•/. 

Icx=  1 : it ; 

Icy=  1 : (it+1) ; 

Jcx=  pml+3 : ( jt+1) -pml-2 ; 

Jcy=  pml+3: jt-pml-2; 

Kc=  pb+1 ; 

field. Ex(Icx, Jcx,Kc)=  0;  '/.Concucting  strip 

field. Ey(Icy, Jcy,Kc)=  0;  '/, Concucting  strip 


Igx=  1 : it ; 

Igy=  1:  (it+1) ; 

Jgx=  1 : (jt+1) ; 

Jgy=  1 : jt ; 

Kg=  (kt+1) ; 

field. Ex(Igx, Jgx,Kg)=  0;  '/.Ground  Plate 

field. Ey(Igy, Jgy,Kg)=  0;  '/.Ground  Plate 


'/,  wall  is  no  longer  on  the  pml 
Iwx=  1 : it ; 

Iwz=  1 :  (it+1) ; 

Jw=  (jt+1) -pml-2 ; 

Kwx=  pb+1 : (kt+1) ; 

Kwz=  pb+l:kt; 

field. Ex(Iwx, Jw,Kwx)=  0;  '/.Concucting  wall 

field. Ez(Iwz, Jw,Kwz)=  0;  '/.Concucting  wall 
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AddPECTF.m 


function  field=  addPECTF (field,  it,  jt,  kt ,  pb,  pml) 

7. 

'/  FUNCTION: 

'/  1.  ASSIGN  ZERO  THICKNESS  PEC  FOR  A  TFW  ANTENNA  TO  BY  ZEROING 

'/  TANGENTIAL  E  FIELDS  CORRESPONDING  TO  CONDUCTOR  LOCATIONS 

'/ 

Icx=  1 : it ; 

Icy=  1 : (it+1) ; 

Jcx=  pml+3 : ( jt+1) -pml-2 ; 

Jcy=  pml+3: jt-pml-2; 

Kc=  pb+1 ; 

field. Ex(Icx, Jcx,Kc)=  0;  '/Concucting  strip 

field. Ey(Icy, Jcy,Kc)=  0;  '/Concucting  strip 

Igx=  1 : it ; 

Igy=  1: (it+1) ; 

Jgx=  1 :  (jt+1) ; 

Jgy=  1 : jt ; 

Kg=  (kt+1) ; 

field. Ex(Igx, Jgx,Kg)=  0;  '/Ground  Plate 

field. Ey(Igy, Jgy,Kg)=  0;  '/Ground  Plate 

Iwx=  1 : it ; 

Iwz=  1 :  (it+1) ; 

Jw=  round(jt/2)+l ; 

Kwx=  pb+1 : (kt+1) ; 

Kwz=  pb+l:kt; 

field. Ex(Iwx, Jw,Kwx)=  0;  '/Concucting  wall 

field. Ez(Iwz, Jw,Kwz)=  0;  '/Concucting  wall 
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AddPECMENZ.m 


function  field=  addPECMENZ (field,  it,  jt,  kt ,  pb,  pml,  dx) 

7. 

7.  FUNCTION: 

7.  1.  ASSIGN  ZERO  THICKNESS  PEC  FOR  A  MENZEL  ANTENNA  TO  BY  ZEROING 

7.  TANGENTIAL  E  FIELDS  CORRESPONDING  TO  CONDUCTOR  LOCATIONS 

7. 


width=  jt-2*pml-4; 
sloty=  round (width/4) ; 
slotx=  round(sloty/2) ; 
pecx=  round(slotx/l)  ; 


7.  WIDTH  OF  STRUCTURE  (CELLS) 

7.  PEC  ON  BOTH  <y>  SIDES  OF  EACH  SLOT 
7.  PEC  BETWEEN  SLOTS  IN  <x> 

7.  SLOT  <x> 


Icx=  1 : it ; 

Icy=  1 : (it+1) ; 

Jcx=  pml+3 :pml+3+sloty ; 

Jcy=  pml+3 :pml+2+sloty; 

Kc=  pb+1; 

field. Ex(Icx, Jcx,Kc)=  0; 

field. Ey(Icy, Jcy,Kc)=  0;  '/Top  long 


Icx=  1 : it ; 

Icy=  1 : (it+1) ; 

Jcx=  (jt+l)-pml-2-sloty: ( jt+1) -pml -2 ; 
Jcy=  jt-pml-l-sloty : jt-pml-2 ; 
field. Ex(Icx, Jcx,Kc)=  0 

field. Ey(Icy, Jcy,Kc)=  0;  /(Bottom  long 


Icx=  1 :pml+2*pecx; 

Icy=  1 :pml+2*pecx+l ; 

Jcx=  pml+3+sloty+l : (jt+l)-pml-sloty-3; 
Jcy=  pml+2+sloty+l : jt-pml-sloty-2 ; 
field. Ex(Icx, Jcx,Kc)=  0; 
field. Ey(Icy,  Jcy,Kc)=  0;  '/First 


Icx=  pml+slotx+2*pecx+l :pml+slotx+3*pecx; 
Icy=  pml+slotx+2*pecx+l :pml+slotx+3*pecx+l ; 
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"/Second 


field. Ex(Icx, Jcx,Kc)=  0; 
field. Ey(Icy, Jcy,Kc)=  0; 

for  n=  1:5 
Icx=  Icx+slotx+pecx; 

Icy=  Icy+slotx+pecx; 
field. Ex(Icx, Jcx,Kc)=  0; 
field. Ey(Icy, Jcy,Kc)=  0;  "/Third  -  Seventh 

end; 

Icx=  Icx(l)+slotx+pecx : it ; 

Icy=  Icy (l)+slotx+pecx : it+1 ; 

field. Ex(Icx, Jcx,Kc)=  0; 

field. Ey(Icy,  Jcy,Kc)=  0;  '/Eighth 

Igx=  1 : it ; 

Igy=  1: (it+1) ; 

Jgx=  1 : (jt+1) ; 

Jgy=  1 : jt ; 

Kg=  (kt+1) ; 

field. Ex(Igx, Jgx,Kg)=  0;  "/Ground  Plate 

field. Ey(Igy, Jgy,Kg)=  0;  "/Ground  Plate 
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ViewData.m 


'/,  function  [A,B]=  ViewData(view) 

7. 

•/.  FUNCTION: 

•/.  1.  READ  DATA  FROM  THW,  TFW,  OR  MENZEL  SIMULATION 

7.  2.  PROCESS  DATA 

•/.  3.  COMPUTE  \ alpha  AND  \beta 

•/.  4.  PLOT  \ alpha  AND  \beta 

7.  5.  SAVE  \ alpha  AND  \beta  TO  EXCEL  SPREADSHEET 

•/. 

clear  all 
close  all 
clc 


•/.  ALLOW  USER  TO  CHOOSE  TRIAL  TO  EVALUATE 
[view,  pathname] =  uigetfile; 
load (  [pathname ,  view]  )  ; 

•/.  DISPALY  THE  RAW  DATA 
for  n=  1 : floor (periods) 
figure 

EZ=  EZDATACn, :) ; 

plot (1 : length (EZ) ,EZ, ’ linewidth’ ,2) 
ylabel(’  E_z  ’); 

xlabel(’  Cell  Number  (x-direction)  ’); 
title([’  Period  #:  ’ ,num2str(n)] ) ; 

grid  on 

end 


•/.  Ez:  GET  ABOUT  TWO  GOOD  WAVES.  IT  MAY  NOT  BE  POSSIBLE  TO  GET  TWO  WAVELENGTHS. 

•/.  ADJUST  ’start’  AND  ’stop’  TO  GET  AS  MUCH  DATA  AS  POSSIBLE.  AT  LEAST  A  HALF 

•/.  WAVELENGTH  IS  REQUIRED  TO  CONTINUE. 

if  f/le9  <6.1  start=  is+40; 

elseif  f/le9  <  7.1  start=  is+10; 

else  start=  is+3;  end; 


if  f/le9  <  6.3  stop=  500; 
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elseif  f/le9  <  6.6  stop=  300; 
elseif  f/le9  <  7.0  stop=  150; 
else  stop=  100;  end; 


•/.  NORMALIZE 

Max=  f ind(abs (Ez (start : start+stop) )==max(abs (Ez (start : start+stop) ) ) ) ; 

Ez=  Ez( (start+Max) : stop+start+Max) /Ez( (start+Max) ) ; 

•/.  Es:  INTERPOLATE  MORE  POINTS  IF  NECCESSARY 
points=  1; 

x=  (1 : length(Ez) )+start+Max-l ; 

xs=  (1 : (1/points) : length (Ez) )+start+Max-l ; 

Es  =  spline(x,Ez,xs) ; 
dxs=  dx/points; 

•/.  USE  LOG(Es)  TO  DETERMINE  ’A’  AND  ’B’ 

logEs=  real (log(Es) ) ; 

figure 

plot (xs*dxs , logEs , ’b’ , ’ linewidth’ ,3) ; 
grid  on 

xlabel ( ’Distnace  from  source  (m) ’ ) ; 
ylabel(’ln  {  E_z  }  (dB)  ’); 

•/.  FIND  NULL  LOCATIONS:  ’nullEz’,  PEAK  LOCATIONS:  ’slopeX’,  AND  PEAK  VALUES:  ’slopeY’ 
slopeX=  1;  slopeY=  0;  nullEs=  0;  down=  1; 
for  nn=  4 : length(logEs) 

if  down  &  logEs (nn) >logEs (nn-1)  &  logEs (nn-1) >logEs (nn-2) 
down=  0 ; 

nullEs=  [nullEs  (nn-2)*dxs]; 

elseif  “down  &  logEs (nn)<logEs (nn-1)  &  logEs(nn-l)<logEs(nn-2) 
down=  1 ; 

slopeX=  [slopeX  (nn-2)]; 
slopeY=  [slopeY  logEs (nn-2) ] ; 

end; 

end 

Lbeta=  nullEs (4) -nullEs (2) ; 

B=  2*pi/Lbeta; 
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peakl=  1; 
peak2=  3; 

A=  (  slopeY(peakl) -slopeY(peak2)  )/(  (slopeX(peak2)-slopeX(peakl))*dxs  ); 

•/.  NORMALIZED  EXPONENTIAL  CURVE  FITTING  ’A’  AND  ’B’ 

EsMatch=  exp(-A*xs*dxs-j*B*xs*dxs) /max(exp(-A*xs*dxs-j  *B*xs*dxs) ) ; 

•/.  COMPARE  DATA  TO  MATCHED  CURVE 

figure 

hold  on 

plot (x,log(Es) , ’r ’ , ’ linewidth’ ,4)  ; 
plot (x, log(EsMatch) ,  ’k’  ,  ’ linewidth’ , 2) ; 
hold  off 
ylabel  (  ’Ez  ’  )  ; 

xlabel ( ’X-direction  cell  number’); 
title(view) ; 

legend ( ’FDTD’ ,  [’ \alpha/k_0=  ’ ,num2str (round (A/k* 100) / 100) , ’ \beta/k_0=  ’ ,num2str (B/k)] ) 
grid  on 

•/.  COMPARE  LOG  DATA  TO  LOG  MATCHED  CURVE 

figure 

hold  on 

plotCx,  (Es) , ’r’ , ’linewidth’ ,4) ; 
plot (x, (EsMatch) , ’k’ , ’linewidth’ ,2) ; 
hold  off 
ylabel  (  ’Ez  ’  )  ; 

xlabel ( ’X-direction  cell  number’); 
title(view) ; 

legend ( ’FDTD ’ ,  [’ \alpha/k_0=  ’ ,num2str (round (A/k* 100) / 100) , ’ \beta/k_0=  ’ ,num2str (B/k)] ) 
grid  on 

•/.  SAVE  DATA  TO  EXCEL  SPREADSHEET  TO  BE  RECORDED 
format  long 

data=  [  f/le9  A/k  B/k  ]; 

N  =  xlsread( ’C : \Documents  and  Settings\Ashley\Desktop\NewSimData’ ,name) ; 

N=  [N ;  data]  ; 

xlswrite( ’C : \Documents  and  Settings\Ashley\Desktop\NewSimData’ ,N,name , ’ A2’ ) ; 
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ViewDataCurve.m 


'/,  function  [A,B]=  ViewDataCurve(view) 

7. 

•/.  FUNCTION: 

•/.  1.  READ  DATA  FROM  CURVE  SIMULATION 

7.  2.  PROCESS  DATA 

•/.  3.  COMPUTE  \ alpha  AND  \beta 

•/.  4.  PLOT  \ alpha  AND  \beta 

•/.  5.  SAVE  \ alpha  AND  \beta  TO  EXCEL  SPREADSHEET 

•/. 

clear  all 
close  all 
clc 


•/.  ALLOW  USER  TO  CHOOSE  TRIAL  TO  EVALUATE 
[view,  pathname] =  uigetfile; 
load (  [pathname ,  view]  )  ; 

•/.  DISPALY  THE  RAW  DATA 
for  n=  1 : floor (periods) 
figure 

EZ=  EZDATACn, :) ; 

plot (1 : length (EZ) ,EZ, ’ linewidth’ ,2) 
ylabel(’  E_z  ’); 

xlabel(’  Cell  Number  (x-direction)  ’); 
title([’  Period  #:  ’ ,num2str(n)] ) ; 

grid  on 

end 

'/,  figure 

'/,  plot  (1 :  length  (angles)  ,  angles ,  ’b’  ,  ’  linewidth’  ,2) 
'/,  title  (’raw  angles’) 

'/,  grid  on 

'/,  xlabel(’cell  number’); 

'/,  ylabel( ’Angle  (degrees)  ’); 

•/.  UNWRAP  THE  ANGLE 
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R=  r+Y-2*dy;  "/,  radius  of  "circle" 

P=  2*pi*R;  7.  perimeter  of  "circle" 

E=  Ez (floor (periods) ; 
clear  Ez 

angles=  unwrap (angles) ; 

'/,  figure 

'/,  plot  (1 :  length  (angles)  ,  angles ,  ’b’  ,  ’  linewidth’  ,2) 

'/,  title ( ’unwrapped  angles’) 

'/,  grid  on 

'/,  xlabel(’cell  number’); 

'/,  ylabel( ’Angle  (degrees)  ’); 

•/.  ORDER  THE  ANGLES 
A=0 ;  Ez=0 ; 

for  aa=  1 : length (angles) 
for  a=  l:length(A) 

if  angles (aa)<  A (a) 

A=  [A(l:(a-1))  angles(aa)  A(a:end)]; 

Ez=  [Ez(l:(a-1))  E(aa)  Ez(a:end)]; 
break; 

elseif  a==length(A) 

A=  [A  angles(aa)]; 

Ez=  [Ez  E(aa)]; 

end 

end 

end 

A(l)  =  []  ;  Ez (1)  =  []  ; 

'/,  figure 

'/,  plot  (1 :  length  (A)  ,  A,  ’b  ’  ,  ’  linewidth’  ,2) 

'/,  title  (’ordered  angles’) 

'/,  grid  on 

'/,  xlabel  (  ’  cell  number  ’  )  ; 

'/.  ylabel( ’Angle  (degrees)  ’); 

•/.  CHANGE  ANGLE  TO  POSITION 

dist=  (A)/360*P;  °/0  vector  of  positions  of  each  point  along  antenna  "arc" 

'/.  figure 
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'/,  plot  (1 :  length (di st)  , dist ,  ’b’  ,  ’  linewidth’  ,  2) 

'/,  title (’ angles  turned  to  distances’) 

'/,  grid  on 

'/,  xlabel(’cell  number’); 

'/,  ylabel(’Arc  distance  to  source  (m)  ’); 

•/.  RAW  LOG  DATA  TO  BE  ANALYZED 
figure 

plot (dist , log(Ez) , ’b’ , ’ linewidth’ ,2) 
title ( ’raw’ ) 
grid  on 

xlabel(’cell  number’); 
ylabel ( ’ I  E_z  I  ’ ) ; 

•/.  FIND  SOURCE 

is=  f ind(log(Ez)==max(log(Ez) ) ) ; 

•/.  Ez:  GET  ABOUT  TWO  GOOD  WAVES.  IT  MAY  NOT  BE  POSSIBLE  TO 

•/.  GET  TWO  WAVELENGTHS.  ADJUST  ’start’  AND  ’stop’  TO  GET  AS  MUCH  DATA  AS 

•/.  POSSIBLE.  AT  LEAST  A  HALF  WAVELENGTH  IS  REQUIRED  TO  CONTINUE. 

if  f/le9  <6.1  start=  is+50; 

elseif  f/le9  <  7.1  start=  is+50; 

else  start=  is+30;  end; 

if  f/le9  <  6.3  stop=  300; 

elseif  f/le9  <  6.6  stop=  300; 

elseif  f/le9  <  7.0  stop=  150; 

else  stop=  200;  end; 

•/.  NORMALIZE 

Max=  f ind(abs (Ez (start : start+stop) )==max(abs (Ez (start : start+stop) ) ) ) ; 

Ez=  Ez( (start+Max) : stop+start+Max) /Ez( (start+Max) ) ; 
dist=  dist (start+Max : stop+start+Max) ; 

'/,  figure 

'/,  plot  (1 :  length  (Ez)  ,Ez,  ’b’  ,  ’linewidth’  ,2) 

'/,  title  (’short  and  normal’) 

'/,  grid  on 
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'/,  xlabel(’cell  number’); 

'/,  ylabel  (  ’  I  E_z  I  ’  )  ; 

•/.  USE  LOG(Es)  TO  DETERMINE  ’A’  AND  ’B’ 

logEz=  real (log(Ez) ) ; 

figure 

plot (dist , logEz , ’b’ , ’ linewidth’ ,2) 
title( ’ log’ ) 
grid  on 

xlabel(’cell  number’); 
ylabel (’In  I  E_z  I  (dB)  ’); 

7.  FIND  NULL  LOCATIONS  ’nullEz’,  PEAK  LOCATIONS  ’slopeX’,  AND,  PEAK  VALUES 
'/,  ’slopeY’ 

slopeX=  dist(l);  slopeY=  0;  nullEz=  dist(l);  down=  1; 
for  nn=  4 : length (logEz) 

if  down  &  logEz(nn)>logEz(nn-l)  &  logEz (nn-l)>logEz(nn-2) 
down=  0 ; 

nullEz=  [nullEz  dist(nn-2)]; 

elseif  “down  &  logEz (nn)<logEz(nn-l)  &  logEz(nn-l)<logEz(nn-2) 
down=  1 ; 

slopeX=  [slopeX  dist(nn-2)]; 
slopeY=  [slopeY  logEz (nn-2)] ; 

end; 

end 

Lbeta=  (nullEz (4) -nullEz (2) ) ; 

B=  2*pi/Lbeta; 

A=  (  slopeY(l) -slopeY(3)  )/(  (slopeX(3) -slopeX(l) )  ); 

•/.  NORMALIZED  EXPONENTIAL  CURVE  FITTING  ’A’  AND  ’B’ 

EzMatch=  exp(-A*dist- j*B*dist) /max(exp(-A*dist- j*B*dist) ) ; 

•/.  COMPARE  DATA  TO  MATCHED  CURVE 

figure 

hold  on 

plot (dist ,Ez , ’r ’ , ’ linewidth’ ,4) ; 
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plot (dist ,real (EzMatch) ,  ’k’ , ’ linewidth’ ,  2)  ; 

hold  off 

ylim(  [-1  1.5]); 

ylabel ( ’Ez ’ ) ; 

xlabel ( ’X-direction  cell  number’); 
title(view) ; 

legend ( ’FDTD’ ,  [’ \alpha/k_0=  ’ ,num2str (round (A/k* 100) / 100) , ’ \beta/k_0=  ’ ,num2str (B/k)] ) 
grid  on 

7.  COMPARE  LOG  DATA  TO  LOG  MATCHED  CURVE 

figure 

hold  on 

plot (dist , log(Ez) , ’r ’ , ’ linewidth’ ,4) ; 

plot (dist , log (real (EzMatch) ) , ’k’ , ’ linewidth’ ,2) ; 

hold  off 

ylabel  (  ’Ez  ’  )  ; 

xlabel ( ’X-direction  cell  number’); 
title(view) ; 

legend ( ’FDTD’ ,  [’ \alpha/k_0=  ’ ,num2str (round (A/k* 100) / 100) , ’ \beta/k_0=  ’ ,num2str (B/k)] ) 
grid  on 

•/.  SAVE  DATA  TO  EXCEL  SPREADSHEET  TO  BE  RECORDED 
data=  [  f/le9  A/k  B/k  time]; 

N  =  xlsread( ’C : \Documents  and  Settings\Ashley\Desktop\NewSimData’ ,name) ; 

N=  [N ;  data]  ; 

xlswrite( ’C : \Documents  and  Settings\Ashley\Desktop\NewSimData’ ,N,name, ’A2’); 


133 


LineSource.m 


function  LineSource (thl , th2 , f ) 

7. 

7.  OUTPUT:  PLOTS  SHOWING  LINE  SOURCE  APPROXIMATION  OF  THE  FAR-FIELD 
7.  PATTERN  OF  THIELE  HALF  WIDTH  ANTENNA  GENERATED  FROM  ALPHA  AND  BETA 

I 

7.  INPUTS : 

7.  \ alpha  AND  \beta  FROM  MEASUREMENTS,  SIMULATION,  AND/OR  TRANS  RES  APPROX 

7.  ’thl’  IS  THE  REFLECTION  ANGLE  FROM  THE  FAR  END 

7.  ’  th2’  IS  THE  REFLECTION  ANGLE  FROMTHE  NEAR  (SOURCE)  END 

7.  ’f’  IS  THE  EXCITATION  FREQUENCY 

7. 

7. 


c=  299792458; 
omega=  2*pi*f ; 
lambda=  c/f; 
kO=  2*pi/lambda; 
uO=  pi*4e-7 ; 
eO=  c~-2/u0; 


7.  speed  of  light  (m/s) 

7.  angular  frequency  (rad/s) 

7,  wavelength  (m) 

7.  wavenumber  (1/m) 

7.  permeability  of  free  space  (H/m) 
7.  permittivity  of  free  space  (F/m) 


load  TRh787wl5e233  7.  GET  TRANSVERSE  RESONANCE  DATA  FOR  ’f’ 

alphaT=  round (AA(f ind (round (f reqs* 100)/ 100==f /le9) ) *1000) / 1000 ; 
betaT=  round (BB (f ind (round (f reqs* 100) / 100==f /le9) ) *100) /100 ; 


if  f==  6 . 7e9 

load  F67  7.  MEASURED  DATA 

az=  F67 ( : ,l)-90; 

mag=  f liplr(F67(  :  ,2) ’ )  ’  ; 

ph=  F67(:  ,3); 

alphaF=  0.0404245257377625;  7.PUT  FDTD  DATA  HERE 
betaF=  0.676644325256347;  7.PUT  FDTD  DATA  HERE 


load  gmz67  7.  MEASURED  DATA 

GMZaz=  gmz67 ( : , 1) -90; 

GMZmag=  (gmz67 ( : , 2) ’ ) ’ ; 

GMZph=  gmz67 ( : ,3) ; 
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elseif  f==  7 . 2e9 


load  F72  */.  MEASURED  DATA 

az=  F72( : ,l)-90; 

mag=  fliplr(F72(: ,2)’)’; 

ph=  F72( : ,3) ; 

alphaF=  0.0251624081283808;  '/.PUT  FDTD  DATA  HERE 

betaF=  0.816086113452911;  '/.PUT  FDTD  DATA  HERE 

load  gmz72  */.  MEASURED  DATA 

GMZaz=  gmz72( : , 1) -90; 

GMZmag=  (gmz72( : ,2) ’ ) ’ ; 

GMZph=  gmz72(  :  ,3) ; 

else 

error (’NO  DATA  FOR  THAT  FREQUENY’); 

end; 

phi=  (0:  .1:180) *pi/180; 

percent=  4/100;  7,  ADD  4°/„  TO  gamma  TO  ACCOUNT  FOR  "LINE  SOURCE/ANTENNA  LENGTH  ISSUE" 
r=  100;  */.  distance  to  far  field  —  2.0  (m)  is  adequate) 

L=  .19;  */.  length  of  antenna  (m) 

•/.  TRANSVERSE  RESONANCE  FORWARD  WAVE 

gammaT=  (betaT* ( 1+percent ) ) *k0- j * (alphaT* ( 1+percent ) ) *k0 ; 

OmT=  L/2* (kO*cos (phi) -gammaT) ; 

ET=  (sin(phi)  .*sin(0mT)  ,/OmT)  ;  */.  NORMALIZED 

maxEtr=  phi  (find  (abs(ET)==max  (abs(ET)  )))  ;  */.  LOCATION  OF  MAIN  BEAM 

ErT=  exp(-alphaT*kO*L)*abs(ET)*exp(-j*thl)  ;  7,  RETURN  WAVE 

Er2T=  exp(-alphaT*kO*L) *ErT*exp(-j *th2)  ;  7.  RETURN-RETURN  WAVE 

Er3T=  exp(-alphaT*kO*L)  *Er2T*exp(-j*thl)  ;  '/.  RETURN-RETURN-RETURN  WAVE 

EtotalT=  abs(ET)  +fliplr(ErT)  +Er2T  tfliplr (Er3T)  ;  7.  ALL  TRANSVERSE  RESONANCE  WAVES 

•/.  FDTD  FORWARD  WAVE 

gammaF=  (betaF* ( 1+percent) ) *kO- j * (alphaF* ( 1+percent) ) *kO ; 

OmF=  L/2* (kO*cos (phi) -gammaF) ; 
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EF=  (sin(phi) . *sin(OmF) . /OmF) ; 

maxEf dtd=  phi (f ind(abs (EF)==max(abs (EF) ) ) ) ; 

ErF=  exp(-alphaF*kO*L) *abs (EF) *exp(-j*thl) ; 

Er2F=  exp(-alphaF*kO*L) *ErF*exp(-j*th2) ; 

Er3F=  exp(-alphaF*kO*L) *Er2F*exp(- j*thl) ; 

EtotalF=  abs (EF)  +fliplr(ErF)  +Er2F  +fliplr (Er3F) ; 


•/.  NORMALIZED 

•/.  LOCATION  OF  MAIN  BEAM 

•/.  RETURN  WAVE 

•/.  RETURN-RETURN  WAVE 

•/.  RETURN-RETURN-RETURN  WAVE 

•/.  ALL  FDTD  WAVES 


•/.  DISPLAY  MEASURED  DATA  AND  FDTD/TR  LINE  SOURCE  APPROXIMATIONS 

figure 

hold  on; 

plot ( (phi*180/pi)+ . 5 ,  20*logl0(EtotalF/max(EtotalF) ) ,  ’r — ’, ’linewidth’ ,  4); 

plot (az , mag-max (mag) , ’k’ , ’ linewidth’ ,3) ; 

plot (GMZaz+2 . 5 ,GMZmag-max(GMZmag) , ’b ’ , ’ linewidth’ ,3) ; 

hold  off; 

grid  on; 

set (gcf , ’Color’ , [1 ; 1 ; 1] ) ; 
set (gca, ’ XDir ’ , ’reverse ’ ) 
axis  (  [0  180  -30  7]); 
title( [num2str(f/le9) , ’  GHz’]); 

legend(’FDTD  /  Line  Source’ , ’Measured  (vias) ’Measured  (tape)’); 

'/,  legend(’Line  Source ’, ’Measured  (tape)’); 
xlabel([’  Angle  from  endfire  (degrees)’]); 
ylabel(’  Normalized  I  E_{H-pol}  I  (dB)’); 

•/.  DISPLAY  BOTH  POLAR  AND  RECTANGULAR 
figure 

subplot (2,1,1) 
hold  on; 

plot ( (phi*180/pi) ,  20*logl0 (EtotalF/max(EtotalF) ) ,  ’b’ , ’linewidth’ ,  3); 
plot (maxEf dtd*ones (20,1) ,linspace(-30, 0 , 20) , ’r ’ , ’ linewidth’ ,  3) ; 
hold  off; 
grid  on; 

set (gcf , ’Color’ , [1 ; 1 ; 1] ) ; 
set (gca, ’ XDir’ , ’reverse ’ ) 
axis (  [0  180  -30  10]); 

title(  [num2str(f/le9) , ’  GHz’,’  \alpha  /k_0=  ’ ,num2str (alphaF) , . . . 
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’  \beta  /k_0=  ’ ,num2str (betaF)] ) ; 

'/,  legend(’FDTD  /  Line  Source  ’,  ’Measured’  )  ; 
xlabel([’  Angle  from  endfire  (degrees)’]); 
ylabel(’  Normalized  I  E_{\theta}  I  (dB)’); 

text (130 ,5 ,  [’  Main  Beam:  ’ ,num2str (maxEfdtd) , ’ ~o  ’ FontSize ’ , 14, . . . 

’BackgroundColor ’ , ’w’ , ’EdgeColor’ , ’r’ , ’Linewidth’ ,3) ; 


subplot (2,1,2) 

p=  f ind(abs (real (20*logl0 (Etotalf )+30) )  <.05); 

polar (phi (p(l) :p(2)) ,  20*logl0(EtotalF (p(l) :p(2) ) /max(EtotalF) )+30) ; 
ylim(  [0  45]  ) 

ylabel ( ’Normalized  I  E_\theta  I  (dB)’); 


•/.  APPROX  MAIN  BEAM 
MAINtr=  acos (betaT) *180/pi 
MAINfdtd=  acos (betaF) *180/pi 

•/.  MAIN  BEAM  DATA 
maxEtr=maxEtr* 180/pi 
maxEf dtd=maxEf dtd* 180/pi 
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TransRes.m 


function  [A,B] =TransRes (f ,h, w,er) 

•/.  WRITTEN  BY:  Dr.  Gary  Thiele 

•/.  MODIFIED  BY:  Greg  Zelinski 

•/.  LAST  MODIFIED:  21N0V2004 

7. 

'/,  INPUT:  ‘f’  is  frequency  (Hz) 

'/,  ‘h’  is  the  substrate  height  (m) 

'/,  ‘w’  is  the  conductor  width  (m) 

'/,  ‘er’  is  the  relative  permittivity  of  the  substrate 

7. 

'/,  OUTPUTS:  ‘A’  is  the  attenuation  constant  normalized  with  free  space  wavenumber 

'/,  ‘B’  is  the  phase  constant  normalized  with  free  space  wavenumber 

7. 

'/,  SOURCE:  K.  S.  Lee,  ‘  'Microstrip  Line  Leaky  Antenna’’,  Ph.D.  diss., 

'/,  Polytechnic  Institute  1986 

c=  299792458;  "/.  speed  of  light 

u0=  pi*4e-7 ;  7,  free  space  permeability 

e0=  c~-2/u0;  70  free  space  permittivity 

k0=  2*pi*f/c;  70  free  space  wavenumber 

k=  kO*sqrt(er);  70  wavenumber  in  medium 

m=  1:50;  70  summation  truncation 

delta=  (er-1) / (er+1) ; 

Q=  sum(  (-delta) . ~m. *log(m)  ); 

syms  kxe  70  symbolic  variable 

kz=  sqrt(k~2  -  kxe~2) ; 

alpha  =  kz/kO ; 

kx=  sqrt (l-alpha~2) *k0 ; 

gamma  =  0.5772-1; 

fe  =  -2*kxe*h/pi* (  (  log(j*kx*h)  tgamma  )/er  +2*Q  -log(2*pi)  ); 

del  =  kz*h/pi*(  (l-er)/er*(  log(j*kx*h)  tgamma  )  +  2*Q  ); 
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chi  =  2*atan(  kz/kxe*tanh(del)  )  -  fe; 


trans=  chi  -kxe*w  +pi ;  */.  TRANSVERSE  RESONANCE  EQUATION  FOR  EH_1 

solkxe=  solve(trans)  ;  "/,  solution  of  symbolic  variable  kxe 

prop=  sqrt(k~2  -  solkxe~2) ; 

B  =  real(prop/kO) ; 

A  =  abs (imag(prop/kO) ) ; 
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MakeTransResFile.m 


'/. 

'/. 

'/. 

'/. 

7. 

'/. 

7. 

•/. 

•/. 

'/. 

7. 

7. 

•/. 


MakeTransResFile  .m 


OUTPUT : 

1.  SAVES  .mat  FILE  TO  THE  CURRECT  DIRECTORY. 

‘freqs’  IS  A  VECTOR  OF  THE  FREQUENCIES  (Hz) 

‘AA’  IS  A  VECTOR  OF  \alpha/kO 
‘BB’  IS  A  VECTOR  OF  \beta/kO 

2.  PLOT  OF  DATA 

NAMING  CONVENTION: 

‘TRh787wl50e233 .mat ’  IS  FOR  AN  ANTENNA  OF 

HEIGHT  787e-6  (m) ;  HALF  WIDTH  7.5  (mm);  AND  DIELECTRIC  CONSTANT  2.33 


clear  all 


close  all 
clc 


format  compact 
tic 

getData  =  0;  '/,  use  "1"  to  get  data  or  "0"  to  just  plot 


if  getData 

range=  (6 :  .  1 : 8 . 2)  *le9 ;  '/.FREQUENCIES 


"/.for  MENZEL’S  antenna: 
h=  . 787*le-3 ; 
w=  15e-3; 
er  =  2.33; 


for  r=  1 : length (range) 
last=  toe; 

[A(r) ,B(r)] =TransRes (range (r) ,h,w, er) ; 

fprintf  (  ’  \nFrequency  '/,.0f  of  '/..Of  took  '/,.lf  s  .  ’  ,r ,  length(range)  ,  toc-last) 

end 


'/.  interpolate  Transverse  Resonance  data  to  1  MHz  increments 
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freqs=  range (1) : le6 : range (end) ; 

AA  =  pchip(range,A,freqs)  ;  "/,  pchip  or  spline 

BB  =  pchip (range, B,f reqs)  ■ 


f ile=  [’TRh’ ,num2str (round(w*le6) )  ,  ’w’ ,num2str (round(w*le4) ) , . . . 

’e’ ,num2str(round(er*100))] 


save(f ile) ; 


else 

f ile=  ’TRh787wl50e233’ ; 
load(f ile) ; 

end 

figure 

plot (f reqs , AA, ’ c ’ , f reqs ,BB , ’r ’ , ’ linewidth’ , 3) ; 

xlabel  (’Frequency  (GHz))’) 

legend (’ \alpha/k_0  ’,’\beta/k_0  ’) 

title  (  [file ,  ’  ’  ]  ) 

grid  on 

set (gcf , ’Color’ ,  [1 ; 1 ; 1] ) ; 

time=  toc/60; 

fprintf  ( ’\n\n\tTransRes  took  2f  minutes  .  \n\n’  ,  time); 
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DataReader.m 


'/,  DataReader.m 

7. 

•/.  PLOTS  A  COMPARISON  BETWEEN  TRANSVERSE  RESONANCE  AND  ONE  OR  MORE  DATA  SETS 

•/. 

clear  all 
close  all 
clc 

load  Compare;  */.  GET  TRANSVERSE  RESONANCE  DATA 

name=  ’180curveS’;  '/,  ENTER  A  DATA  SET  TO  DISPLAY 

N  =  xlsread( ’C : \Documents  and  Settings\Ashley\Desktop\NewSimData’ ,name) ; 
f =  N( : , 1) ;  a=N(:,3);  b=N(:,6); 

name2=  ’Single’;  '/  ENTER  ANOTHER  DATA  SET  TO  DISPLAY 

N2  =  xlsread( ’C : \Documents  and  Settings\Ashley\Desktop\NewSimData’ ,name) ; 

f2=  N2( : , 1) ;  a2=N2(:,2);  b2=N2(:,6); 


figure 
hold  on 

plot (f reqs , AA, ’r ’ , ’ linewidth’ ,4) ; 
plot (f 2 , a2 , ’k’ , ’ linewidth’ ,2) ; 
plot (f , a, ’b . ’ , ’Marker Size ’ ,20) ; 
plot (f reqs ,BB , ’r ’ , ’ linewidth’ ,4) ; 
plot (f 2 ,b2 , ’k’ , ’ linewidth’ ,2) ; 
plot (f ,b, ’b . ’ , ’Marker Size ’ ,20) ; 
hold  off 


axis ([5.6  8.2  0  1.1]); 
set (gcf , ’Color’ ,  [1 ; 1 ; 1] )  ; 
xlabel ( ’Frequency  (GHz)’); 

text (6 . 7 , . 85 , ’ \beta  /k_0’ , ’BackgroundColor ’ , ’ w’ , ’FontSize’ , 16) ; 
text (7, . 15 , ’ \alpha  /k_0 ’ , ’BackgroundColor’ , ’w’ , ’FontSize ’ ,16) ; 
legend ( ’Transverse  Resonance ’,’ Straight  (single)’,... 

’180~o  Curve  (single) ’, ’Location’ , ’East ’) ; 


grid  on 
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was  demonstrated  that  shorten  production  time  and  improved  the  far-field  pattern. 
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